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UNSTEADY HOVERING WAKE PARAMETERS IDENTIFIED 
FROM DYNAMIC MODEL TESTS 

Part I of Final Report under Contract NA£2*7613 

Abstract 

This paper reports on the development of a 4-bladed model rotor 
that can be excited with a simple eccentric mechanism in progressing 
and regressing modes with either harmonic or transient Inputs. 
Parameter Identification methods were applied to the problem of 
extracting parameters for linear perturbation models. Including rotor 
dynamic Inflow effects, from the measured blade flapping responses to 
transient pltcb-stlrrlng excitations. These perturbation models were 
then used to predict blade flapping response to other pitch-stirring 
transient inputs, and rotor wake and blade- flapping responses to 
harmonic inputs. The viability and utility of using parameter 
identification methods for extracting the perturbation models from 
transients are demonstrated through these combined analytical and 
experimental studies. 



Preface to Final Report under Contract NAS2-7613 


Work under Contract MS2-7613 started on July 1, 1973. The contract 
was originally awarded for a 3 year period. 

D(« to the slower than anticipated progress of the experimental work, 
not all research goals had been achieved by 30 June 1976. Since less than 
the anticipated cost for personnel and equipment had been spent, the 
research contract was extended by a year without Increase In funding. 

The research goals as stated In the contract were: 

(a) Assess analytically the effects of fuselage motions on stability 
and random response. The problem is to develop an adequate but not 
overly complex flight dynamics analytical model and to study the 
effects of structural and electronic feedback, particularly for 
hingeless rotors. 

(b) Study by computer and hardware experiments the feasibility of ade- 
quate perturbation models from non-linear trim .ondltlons. The 
problem is to extract an adequate linear perturbation model for the 
purpose of stability and random motion studies. The extraction la 
to be performed on the basis of transient responses obtained either 
by computed time histories or by model tests. 

(c) Extend the experimental methods to assess rotor wake-blade 
Interactions by using a 4-bladed rotor model with the capability 
of progressing and regressing blade pitch excitation (cyclic pitch 
stirring), by using a 4-bladed rotor model with hub tilt stirring, 
and by testing rotor models in sinusoidal up or side flow. 
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Including the final report^ 10 reports under Contract NAS2-7613 
have been submitted. They are listed as P. 1 to P. 10 at the end of 
the Preface. P. 1 and P. 10 pertain to research goal (a). P. 2, P. 4, 

P. 6, P. 7, P. 8, P. 9, pertain to research goal (b). P. 3 and P. 5 
pertain to research goal (c). The latter is not as yet complete since 
neither hub tilt stirring nor testing is sinusoidal up or side flow 
has been performed. While P. 10 describes only work done during FY 1977, 

P. 8 and P. 9 combine both FY 1977 work results and summaries of earlier 
results, so that the three parts of the Final Report can be read without 
recourse to the earlier reports. P. 8 includes much new material not 
available when the preceding Yearly Report P. 7 was written. The 
experimental data of P. 9 have all been obtained in FY 77. 

So far 3 publications came out of the research under Contract NAS2-7613. 
They are listed as P. 11, P. 12, P. 13. 

List of Reports and Papers 
under Contract NAS2-7613 

P 1. Hohenemser, K. H. and Yin, S. K. , "Method^ Studies Toward Simplified 
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paration of System Identification from Transient Rotor Model 
Teats", Part II of First Yearly Report under Contract NAS2-7613, 

June 1974. 

P 3. Hohenemser, K. H. and Crews, S. T., "Experiments with a Four-Bladed 
Cyclic Pitch Stirring Model Rotor", Part III of First Yearly Report 
under Contract NAS2-7613. 

P 4. Hohenemser, K. H. , Banerjee, D. and Yin, S. K., "Methods Studies 
on System Identification from Transient Rotor Tests", Part I of 
Second Yearly Report under Contract NAS2-7613, June 1975. 



11 


P 5. Hohenamser, K. H. and Crews, S. T., ''Additional Experiments with a 
Four*Bladed Cyclic Pitch Stirring Model Rotor". Part II of 
Second Yearly Report under Contract MAS2-7613, June 1975. 

P 6. Hohenemser, K. H. , Banerjee, D. and Yin, S. K., "Rotor Dynamic 

State and Parameter Identification from Simulated Forward Flight 
Transients", Part I of Third Yearly Report under Contract 
NAS2-7613, June 1976. 

P 7. Hohenemser, K. H. and Crews, S. T., "Rotor Dynamic State and 

Parameter Identification from Hovering Transients", Part II of 
Third Yearly Report under Contract NAS2-7613, June 1976. 

P 8. Hohenemser, K. H. and Crewe, S. T. , "Unsteady Hovering Rotor 

Wake Parameters Identified from Dynamic Model Tests", Part I of 
Final Report under Contract NAS2-7613, June 1977. 

P 9. Hohenemser, K. H. and Banerjee, D. , "Application of System 
Identification to Analytic Rotor Modeling from Simulated and 
Wind Tunnel Dynamic Test Data", Part II of Final Report under 
Contract NAS2-7613, June 1977. 

P 10. Hohenemser, K. H. and Yin, S. K. , "Finite Element Stability 
Analysis for Coupled Rotor and Support Systems", Part III of 
Final Report under Contract NAS2-7613, June 1977. 

P 11. Hohenemser, K. H. and Yin, S. K. , "On the Use of First Order 

Rotor Dynamics in Multiblade Coordinates", 30th Annual National 
Forum of the American Helicopter Society, May 1974, Preprint 831. 

P 12. Banerjee, D. and Hohenemser, K. , "Optimum Data Utilization for 
Parameter Identification with Application to Lifting Rotors", 
Journal of Aircraft, Vol. 13, No. 12, December 1976, pp. 1014-1016. 

P 13. Banerjee, D. , Crews, S. T. , Hohenemser, K. H. and Yin, S. K., 
Identification of State Variables and Dynamic Inflow from Rotor 
Model Dynamic Tests", Journal American Helicopter Society, 

Vol. 22, No. 2, April 1977. 



Hi 


TABLE OF CONTENTS 

No. Page 

1. Introduction 1 

1.1 Scope of Work 1 

1.2 Related Work 2 

2. The Four Bladed Rotor Model 5 

2.1 The Rotor Head 5 

2.1.1 Rotor Head Properties 6 

2.1.2 Rotor Head Adjustments 9 

2.2 Rotor Test Set--Up 13 

2.2.1 Test Stand 13 

2.2.2 Analog Measurement Systems and Control Units • 16 

2.3 Calibration and Test Procedure 21 

2.3.1 Harmonic Test Calibration and 

Testing Procedure 22 

2.3.2 Transient Test Calibration and 

Test Procedure 22 

2.4 A Comparison of Test Procedure and Test 

Set Physical Conditions 24 

3. Data Processing 27 

3.1 The Mini-Computer System 27 

3.2 The Data Reduction Algorithms • • . • 28 

3*2.1 Sampling Programs 28 

3.2.2 Calibration Programs 31 

3.2.3 Reduction Programs and Diagnostics 32 

4. Analytical Treatment of the Hovering Unsteady 

Rotor Wake Coupling 36 

5. Development of the Identification Algorithm 42 

5.1 Theoretical Formulation 43 

5.2 Practical Considerations 44 

5.3 Implementation of the Algorithm 47 

5.3.1 Measurements 47 



Iv 

5.3.2 Choice of Parameters 48 

5.3.3 Sensitivity Equations and Matrix SO 

5.3.4 Initial Conditions 51 

5.3.5 Measurement Bias 52 

5.4 Program Description 53 

6. Results and Discussion 55 

6.1 Harmonic Tests 55 

6.2 Data Length Studies 66 

6.3 Identified Parameters 73 

6.4 Predictions 85 

7. Conclusions 92 

7.1 Conclusions Regarding Method . 92 

7.2 Conclusions Regarding Rotor Dynamic Inflow Effects . . 92 

8. References 94 

9. Appendices 98 

9.1 Notation 98 

9.2 First Order Unsteady Wake Identification Program . . . 102 



V 


LIST OF FIGURES 

No. Page 

1. Rotor Head Schematic 6 

2. Rotor Head Photograph 7 

3. First Node Bending Moments 9 

4. Collective Pitch Adjustment 10 

5. Rotor Hovering Test Stand 12 

6. Harmonic Pitch Stirring Excitation System • 14 

7. Transient Pitch Stirring Excitation System 14 

8. Short Transient Pitch Stirring Excitation System ... . 14 

9^ Instrumentation System 18 

10. Hot Wire Anemometer Probe Stand 20 

11. Resolver Normalization 31 

12. Transient Pitch Resolver Sample Window 34 

13. Downwash Response with Harmonic Excitation versus , 38 

14 a-d 0Q =• 5°, = 1.17 Harmonic Test Downwash Measurements . . 58-61 

15. = 5°, = 1.17, 0 ) * .933 Equivalent Linear Dynamic 

Downwash 62 

16. 9^ = 2°, = 1,17 Harmonic Stirring Test 63 

17. 0^ = 5°, oij^ = 1.17 Harmonic Stirring Test ........ 64 

18. 0^ * 8°, = 1.17 Harmonic Stirring Test 65 

19. A' Data Length Study 68 

20. L, T Data Length Study 69 

21. A, L, T Data Length Study 70 

22 a,b P2, A, L, T Data Length Study 71-72 

23 a A' and Fit Factor versus 0 for the Quasisteady 

Wake Model ° 77 

23 b Conq>arison Between Test and Theoretical Results for 

the Quasisteady Wake Model 78 



Vi 

24 a Fit Factor » RR veraua 6 for the Two Parameter 

Unsteady Wake Model . ? 79 

24 b L and T versus 6 for the Two Parameter 

Unsteady Wake Model 80 

25 a Fit Factor and A versus 6 for the Three 

Parameter Unsteady Wake mdel 81 

25 b L and T versus Oq for the Three Parameter 

Unsteady Wake Model 82 

25 c Theoretical Values of L and t versus 8 83 

o 

26. A, L, and t versus for the Unsteady Wake Model .... 84 

27 a,b Progressing Pitch Stirring Transient Test Prediction . . . 87-88 

28. Regressing Test Prediccion Slow Transient 89 

29. Regressing Test Prediction Fast Transient 90 

30. Progressing Test Prediction Short Transient 91 



vll 

LIST OF TABLES 

No. Page 

1. Natural Frequencies at Test Eotor Speeds 9 

2. Calibration and Test Procedure - Transient Tests , 23 

3. Time Comparisons of Test Methods 25 

4. Physical Conditions of Tests 26 

5. Data Reduction 29 

6. Parameter Choice Method Studies 49 

7. Average Inflow at h/R =1.28 56 

8. Averaged Fit Factors . 73 

9. Transient Identification Comparisons 76 



UNSTEADY HOVERING WAKE PARAMETERS IDENTIFIED 
FROM DYNAMIC MODEL TESTS 


1. INTRODUCTION 

While parameter identification from transients is a relatively 
well developed method to extract rotor derivatives from rotorcraft 
flight tests (3ee References 1, 2 , 3, 4, 5, and 6) dynamic rotorcraft 
or model testing has been restricted to frequency response testing 
(References 7 , 8 , 9 , 10, 11, 12, and 13), limit cycle searches 
(Reference 14), and studies of the decay of damped modes (References 
15, 16, and 17). Frequency response tasting implies exciting the 
rotor with a series of discrete constant frequency Inputs over a range 
of frequencies. A rotor transient test that sweeps this frequency 
range takes much less testing time. The question then is whether 
parameter Identification techniques are suitable to extract information 
from transient tests t.iat is equivalent to the information that can 
be extracted from the harmonic tests. 

1.1 SCOPE OF WORK 

The work described herein has two goals. One goal was to build a 
rotor model capable of both harmonic and transient testing, use 
parameter identification methods to extract information from the 
transient tests, use traditional frequency response methods to extract 
Information from the harmonic tests, and compare the two methods of 
testing and data evaluation. Characterization of the hoveling rotor 
dynamic wake Is the second major goal. Some experimental work toward 
this end is presented in Reference 11, which also contains a mathemati- 
cal model for fitting the experimental results. Reference 18 presents 
a theory for determining the parameter values of this mathematical 
model and also extends the model to forward flight. 

A four bladed rotor model was designed to study the structure of 
the rotor dynamic wake. The rotor model can be excited in progressing 
and regressing modes with either harmonic or transient cyclic pitch 
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stirring inputs. Collective pitch can be changed from test to test 
but not during a test run* VIhen harmonic inputs were used, rotor 
quantities measured were pitch input position, rotor shaft position, 
one blade flap^bending deflection, and rotor wake velocities at a 
variety of radial positions. The deflection and velocities were 
Fourier analyzed at multiples of the excitation frequency* The Fourier 
coefficients could be compared to coefficients generated by an analyti** 
cal model* All analytical models are perturbation models about a 
particular trim condition* When transient inputs were used, rotor 
quantities measured were pitch input position, rotor shaft position, and 
four blade flap-bending deflections* Parameter Identification methods 
were used to find the parameter values for analytical models from these 
measurements. The information obtained from both types of experiments 
could be compared by using the transient test generated parameters to 
obtain Fourier coefficients from the analytical models that could then 
be compared to the harmonic test Fourier coefficients. 

Parameter identification methods reduced model testing time by a 
factor of more than 80 for each collective pitch setting studied, and 
produced more easily interpretable results with respect to the analyti- 
cal models studied* 

1.2 RELATED WORK 

Parameter identification techniques have not been applied to 
rot^rcraft model tests before, so there exists no related work in that 
area. However, such techniques have been applied to helicopter flight 
testing (References 1-6). The tests in References 1-4 were used to 
determine derivative models for the maneuvers performed in the tests. 

One use for these models could be In constructing pilot training 
simulators for the particular helicopter that was tested. Another use 
could be to substantiate or improve the analysis, though this use was 
not emphasized in the cited references. In contrast parameter values 
identified from the tests presented in the following are directly 
compared with mathematical model parameter values. In other words, the 
parameter identification techniques are used here to directly explore 
the relationship between mathematical model and the physical dynamic 
system. 

The oldest mathematical rotor models were baaed on the actuator 
disc concepts developed by Glauert and Lock in the 1920's and used 
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loomentusi theory to derive values for the induced inflow through the disc 
and its effects on rotor thrust* These models did not attempt to 
describe the rotor structural motions (rotor dynamics) or transient 
thrust effects, except for the damping from steady pitch or roll rates. 
These models assumed a static uniform wake structure, the effects of 
the dynamic wake on rotor dynamics was Ignored. Steady state models 
and their modem day extensions give fairly good values of gross rotor 
performance. 

Recently much effort has been devoted to the study of rotor wakes, 
primarily for performance prediction. A good survey of the work done 
in the Western world is presented in Reference 19. The Russian treat- 
ment of this same subject appears in Reference 20. These works lean 
heavily toward vortex theory in their treatment of the wake. Vortex 
theories appear to be quite promising for predicting quasi-steady 
rotor wake-structural coupling at moderate to high advance ratios 
(References 20, 21) where certain simplifying assumptions can be made 
about the distribution of wake vorticity. It is still an open question 
as to how useful general vortex theory will become for predicting 
dynamic rotor wake-structural coupling because of the complexity of the 
problem. 

Theoretical and experimental studies of the effects of steady wake 
asymmetries have shown that wake-structural coupling is responsible 
for deficiencies in control effectiveness and for changes in static 
derivatives not only at low but at high advance ratios (References 22, 
23, 24, 25). Reference 24 presents a momentum theory formulation for 
predicting the quasi-steady rotor wake-structural coupling and compares 
this theory to experiment, vortex theory, and an empirical (best fit 
to the data) wake. 

The classical work by Loewy, Reference 26, on the unsteady 
aerodynamics of rotors, represents an extension of the Theodoraan and 
Kussner theory of oscillating airfoils in plane flow to the lifting 
rotor in the substall region at zero advance ratio and predicts an 
appreciable reduction in the aerodynamic damping of flap-bending modes 
in regions where the flap-bending frequency is an Integral multiple of 
the rotation frequency. The Loewy theory was verified in experiments 
conducted by Oaughaday et al, Reference 13. The Loewy effect la 
substantial only at low rotor lift. 
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Carpenter in 1953, Reference 15, Introduced dynaialc Induced down- 
wash in a study of the effect of transient collective pitch changes on 
thrust behavior. Carpenter et al were apparently the first and until 
now only group to make direct measurements of a dynamic wake, although 
many measurements have been made of quasi-steady wakes (see Reference 
19 for a partial literature survey). They showed that for rotors 
subject to collective pitch xavssp inputs, rotor behavior could be 
predicted by including an unsteady wake which is based on the concept 
of a participating air mass. This air mass is equivalent to the air 
mass that participates in the axial acceleration of a circular disk of 
the same radius as the rotor. 

The little that is known to date about cotpled rotor wake and 
cyclic flap-bending characteristics has been deduced from the rotor 
flapping frequency response to cyclic pitch input tests reported In 
References 7-10. Correlations of such test data with a simple 
dynamic rotor wake model based upon a time delayed momentum wake are 
given in Reference 11. Theoretical values of the wake time constants 
are determined in Reference 18 from potential flow around a solid disk. 
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2. THE FOUR BLADED ROTOR MODEL 

The four bladed rotor model was designed specifically for the 
purpose of studying rotor flapping response to dynamic pitch stirring 
inputs. A dynamic pitch stirring input is a cyclic pitch input where 
the position of the maximum pitch input angle either progresses around 
the rotor azimuth angle in the direction of rotor rotation, or 
regresses around the rotor azimuth angle opposite to the direction of 
rotor rotation. The complicated control system needed to control a 
rotor in flight is not needed to produce dynamic pitch stirring inputs. 

In fact a flight control system in a model rotor, which is built to study 
a particular aspect of flight dynamics, is a liability. The many bearings 
necessary to achieve a flight control system usually produce excessive 
dead areas, play, or slop in control systems of model rotors. Our model 
rotor has only six bearings, four pitch feathering bearings, one for 
each blade, and two cyclic pitch stirring input bearings. There is no 
measurable play in the control system. The rotor was designed and built 
to minimize rotor pitch-flap and flap-lag structural coupling over the 
range of experimental test conditions. There is no structural inter- 
blade coupling because the hub is very stiff compared to tne blades. 
Therefore the only possible inter-blade coupling is through the wake 
structure. There are no test stand natural frequencies in range of the 
test excitation frequencies to Interfere with test results. Flapping 
response to cyclic pitch stirring Inputs for this model rotor have been 
effectively isolated from structural coupling with other rotor structural 
motions. 

2.1 THE ROTOR HEAD 

Figure 1 is a schematic of the rotor head and Figure 2 is a photo- 
graph of the rotor head. The blades are driven by the hub, which is 
screwed into the rotor shaft RS, through flanged extended inner ring 
micro-precision bearings. The cyclic pitching motion is generated by 
an inner shaft I whose rotation can be Independent of the rotor shaft 
RS. Relative rotation of the inner shaft I to the rotor shaft RS 
produces via eccentric E, bending flexures B and pitch shafts S 
feathering oscillations of the pair of opposing blades connected with 
flapping flexures F to pitch shafts S. The amplitude of cyclic pitch 
transmitted to the blades is proportional to the eccentricity of the 
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rod in the eccentric E. When the inner shaft is not rotating, a 
steady cyclic pitch excitation is transmitted to the blades. The 
difference in angular velocities of the eccentric shaft 1 and rotor 
shaft RS then produces progressing or regressing cyclic pitch excitations. 
If I rotates in the same direction as RS, then the cyclic pitch is 
progressing in the direction of rotor rotation, if I rotates in the 
opposite direction from RS, then the cyclic pitch is regressing. The 
bearings attached to the bending flexures B will take a plus or minus 
3 ° misalignment. All tests were conducted with a 1.4° eccentric. 

2.1.1 Rotor Head Properties 

The rotor was designed to fit into a wind tunnel with a 24" by 24" 
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cross-section. A rotor diameter of 16.5" was chosen. The NACA 0012 
fiberglass coated styrofoam blades are described In References 27 and 
28. The blade cord Is 1" yielding a rotor solidity ratio of .154. 

Measurements to determine the structural properties of the blade 
and the blade-hub attachment were made with a stationary rotor. The 
blade fixture nut (Figure 2) was excited with a variable frequency 
magnetic field. The blade flapping* edgewise (or lead-lag), and 
torsion (or pitch) response were monitored by appropriate wiring of the 
4 flapping strain gages (Figure 2), and resonant frequencies were 
determined. The measuroients gave 14.0 Hz and 166 Hz for the first and 
second flap bending natural frequencies, 180 Hz and 271 Hz for the first 
and second torsion natural frequencies and 114 Hz for first lag 
bending natural frequency. 

The blade and blade hub attaclunent mass distribution was measured. 

A finite element analysis (described In References 27 and 29) with 20 
elements used the experimentally determined non-rotating first and 
second natural frequencies to construct the unknown flap bending stiff- 
ness distribution of the blade and blade bub attachment . The analysis 
was then used to determine the first and second natural flap bending 
frequencies during rotation, the mode shapes, and the flap bending 
moments. The computed bending moments per unit tip deflection of the 
first mode were found to be independent of rotor velocity at the center 
of the bending flexure strain gages (Figure 3). The value .573 Ibf-ln./ln. 
tip deflection could be used Independent of rotor velocity. This value 
was checked experimentally by applying simulated centrifugal forces. 

For both analysis and experiment the flapping angle 3 Is arbitrarily 
defined by the straight line from the rotor center to the first mode 
deflection at .7 radius. One Inch tip deflection represents 6.32*^ 
flapping angle. The flap bending moment of the strain gage per degree 
of flapping angle is .091 In. -Ibf /degree. 

Tests were conducted at rotor rotational frequencies (n) of 
approximately 21, 27, and 36 Hz (Table 1). The maximum frequency of 
excitation was 1.5 n. Because of the substantial frequency spread 
between excitation and the higher natural frequencies (Table 1), second 
mode flap bending, elastic flap torsion coupling and elastic flap lag 
coupling were Ignored in the dynamic theories used to predict experi- 
mental results. 
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Table 1 

Natural Frequencies at Test Rotor Speeds 


rotor 

1st natural 

2nd natural 1st natural 

1st natural 

maximum pitch 

speed-n 

flap bending flap bending 

torsion 

edgewise 

excitation 

21 Hz 

1.26 n 

8.1 n 

8.6 n 

5.4 n 

1.5 n 

27 Hz 

1.18 n 

6.4 n 

6.7 n 

4.2 n 

1.2 n 

36 Hz 

1.12 n 

4.9 n 

5.0 n 

3.2 n 

1.2 n 

All 

quantities used to calculate 

the blade 

Lock number 

Y are known 


within three digits except for the blade lift slope a« The blade 
Reynolds number at *7 radius varies from .3 to .6 x 10^ through the test 
range. Published test data (References 30 and 31) indicate that 
the 0012 airfoil should have a lift slope of about 5.7 at low angle of 
attack near these Reynolds numbers. However different sources of 
published data on the same airfoil can show characteristics that differ 
by as much as 5 percent (Reference 32). These differences can 
probably be attributed to surface differences between the respective 
airfoils tested. The Lock number for the test blades Is 3.6 t .2. 

2.1.2 Rotor Head Adjustments 

The blades are attached to the hub at the blade root fixtures 
shown in Figure 2. There is a cutout on one side of the blade roo^ 
that allows the blade to slip Into the fixture channel. Two screws are 
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used to bold the blade to the fixture. Blade pre-cone can be eliminated 
by trial and error shimming of the blade root to the fixture nut. Pre- 
lag can also be eliminated by careful blade placement before tightening 
the screws. A special fixture was built to facilitate the elimination 
of blade pre-lag. 

Collective pitch is adjusted by locking the feathering shafts S 
(Figure 1) into different positions with respect to the eccentric 
bending flexures B by tightening the feathering nuts N. Collective 
pitch can only be adjusted with the rotor stopped between test runs. 
Relative collective pitch changes are measured with an optical device. 
Figure 4. 



light 

source 


Figure 4 Collective Pitch Adjustment 

Light strikes the mirror (Figure 2) on the blade flexure root and 
reflects on the measuring scale. With this device collective pitch 
changes on each blade can be made to accuracies of ±.05^. The 
reference value of the collective pitch of each blade, which was used 
for the base of the blade pitch changes, is difficult to find. Four 
sets of hovering tests were conducted. Between each set of tests, 
refinements were made in setting the reference value of collective 
pitch for each blade. Table 4 in Section 2.4 defines the physical 
condition of the rotor for the four sets of tests. All blade pitch 
changes were carefully recorded. 

The rotor shaft is outfitted with pieces of refleculng tape at 
the azimuthal positions of each blade. A photo-electric pickup causes 
a stroboscope to fire as each blade passes the pickup. The viewer sees 
the 4 blade tips (each a different color) at the same vertical special 
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location. Tracking the rotor-getting all the blades to rur in the same 
plane - is achieved by adjusting the individual blade pitch settings 
until all the blade tips are super- Imposed upon each other. The rotor 
can be tracked to an accuracy of greater than .1^ in flapping angle. 

Frequency response tests were conducted first. During these tests » 
downwash velocities were measured with a hot wire anemometer. Keeping 
the rotor tracked, the point of near zero flow velocity was found by 
adjustment of the rotor collective pitch setting. Tests were run at 
plus 2 , 5, and 8 degrees pitch from this point with a pitch stirring 
excitation of 1.4^ and with the rotor 1.28 radii above the ground 
plane. A steady cyclic pitch input moves each blade through one com- 
plete sinusoidal pitching cycle per revolution. Therefore, if the 
perturbations are linear, the blade flapping response should contain 
only Fourier components of once per revolution. During analysis the 
second Fourier coefficient of the flapping response was found to be 
approximately 10 percent of the first. Careful measurement revealed 
that some of the blades had pre-cone angles other than ^ero degrees. 

The spread was approximately .3 degrees. Because of the pre-cone 
spread the tracked blades also had a collective pitch spread of ±.7^. 

The collective pitch given for the frequency response tests is an 
average value. 

The pre-cone was then corrected and transient tests set 1 were 
conducted at collective pitch settings of 1.5, 4.5, and 75 degrees. 
During analysis, the second Fourier coefficients of the blade flapping 
responses to a one per revolution pitch excitation were found to be 
6 to 8 percent of the first. Careful measurements revealed that some 
of the blades had pre-lag angles other than zero. The spread was 
approximately .7^. Because of the geometry of the blade root attach- 
ment, a pre-lag causes a pre-cone. The collective pitch spread for 
this set of tests was ±.4^« 

The pre-lag and its respective pre-cone were then corrected and 
transient tests set 2 were conducted at collective pitch settings of 
3.3, 6.3, and 9.3 degrees. During analysis the second Fourier 
coefficients of the blade flapping responses to a 1 per revolution 
pitching excitation were now found to be 2 to 4 percent of the first 
Fourier coefficients rather than the 6 to 8 percent found in transient 
tests set 1. 






rotor 


^Ilp rin 
assembly 






13 


At this point the rotor was put Into the wind tunnel at zero shaft 
angle of attack and at zero collective pitch setting* This position 
yielded a near zero flapping response from all blades at high advance 
ratio and zero cyclic pitch. Thus the zero collective pitch baseline 
determined previously was confirmed. 

The rotor was then removed from the wind tunnel and transient 
pitch stirring tests set with 5 degree collective pitch angle » were 
run. 

2.2 ROTOR TEST SET-UP 

The rotor test set-up consists of the test stand and the analog 
measurement system. The design of the measurement system was In part 
the result of decisions that had to be made as to how much of the 
signal processing and measurement functions should be done with analog 
devices and how much should be done with digital algorithms on the 
various mini-computers used for data processing. Signal processing and 
measurement functions were Implemented on the digital side whenever 
possible in order to reduce hardware costs. The digital measurement 
functions are described In detail In Section 3. 

2.2.1 Test Stand 

The rotor hovering teat stand can be seen in Figure 5. The ground 
height is adjusted by moving the ground plate support rods (Figure 5). 
The rotor height above the ground plate for each series of tests is 
given In Table 4 section 2.4. The ground plate was constructed to 
achieve C^s far as possible) an azlmuthally symmetric test condition. 
The rotor center Is 4 rotor radii from the edge of the ground plate. 

The ground plate support legs are 2.8 rotor radii away from the rotor 
center and present an obstruction over 5 percent of the circumference. 
The nearest object » apart from the ground plate, to the rotor Is the 
celling which is about 9 radii away from the rotor center. There Is no 
air circulation system in the room. With the windows closed, local air 
velocities are near zero. 

The ro.:or Is leveled by adjusting the lengths of the 4 legs of 
the rotor stand. The rotor is level when a rotor blade has the same 
height above the ground plate (Figure 5) at all azimuth angles. The 
rotor stand Is described In References 27 and 33. The only major 
change that has been made to the rotor stand is the replacement of the 
12 slip ring assembly by a 20 slip ring assembly. 
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The pitch ertlrrlng condition can ba described by 
e - sin (i|»p(t)) 

0JJ “ cos (4ip(t)) ■' 

where *4)p is the inner shaft azlomthal position and is defined so that 
6^ and 6^^ are forward and left cyclic pitch inputs respectively. 

Three cyclic pitch stirring excitation systems can be mounted on the 
rotor stand. 

The harmonic cyclic pitch stirring excitation system is shown in 
Figure 6 and is described extensively in References 27 and 29. The 
drive motor turns the rotor shaft RS (figure 1) through a timing belt 
drive. The rotor shaft drives the inner shaft 1 (Figure 1) through a 
reversible gear train (Figure 6). The gear ratios can be changed. Due 
to the reversing gears, the inner shaft I can be driven in the same 
direction as the rotor shaft, thus providing a progressing pitch 
stirring excitation, or opposite to the direction of rotor shaft rotation, 
thus providing a regressing pitch stirring excitation. Gear ratios are 
available for providing cyclic pitch stirring excitations in a space 
fixed reference system of a ±.067, ±.l, ±.2, ±.4, j:.6, ±.8. 

A transient pitch stirring excitation system that makes use of a 
variant of the double motor drive described in References 27 and 28 is 
shown in Figure 7. The inner shaft (Figures 1 and 7) is driven by a 
variable speed reversible motor through a 1 to 8 speed ratio timing 
belt drive C^^lgure 7). When the inner shaft motor is stationary the 
rotor cyclic pitch is constant with a space fixed excitation frequency 
of 0. Engaging the inner shaft motor accelerates the inner shaft to a 
preselected constant speed and produces either a progressing or regressing 
constant frequency value (depending upon the direction of motor rotation). 
The transient between these two steady values is used for parameter 
Identification. This transient takes the form 

hi ’ %o 

where is a constant. 4'p(t) is zero until the transient starts, the 
inner shaft motor is engaged. *^’p(t) is monotonlcally increasing during 
the transient, as the motor accelerates, 'l'p(t) is constant after the 
transient is over, the motor reaches constant speed. Friction between 
the rotor shaft RS (Figure 1) and the inner shaft I will drive the 
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inner shaft (and Che Inner shaft motor) In the direction of rotor 
rotation. In transient tests sets 1 and 2 (Table 4), the inner shaft 
was restrained by hand until the inner shaft motor was engaged. In 
transient tests set 3, Che inner shaft was released by a solenoid 
(Figure 7) Chat frees the inner shaft for rotation at the same time 
the inner shaft motor is engaged. 

A second method of generating a transient cyclic pitch excitatluu 
is presented in Figure 8. This device has a driver arm attached Co the 
inner shaft which is driven by a spring Chat wraps around a mandrel 
which is below the inner shaft. A solenoid is engaged that releases 
Che driver arm. The spring accelerates Che driver until it hits a 
catcher which is designed to eliminate rebound of the driver arm. The 
catcher can be positioned so that it allows an inner shaft travel of 
90°, 180°, or 270°. The spring torque can be Increased by winding it 
up more. This excitation takes Che form 
- sin («|; + i|< (t)) 

011 = cos C«l<p^ + 1 <p(t)) 

where ip is a constant. iJ; (t) » 0 before the transient starts, before 
po p 

the solenoid releases the driver arm. ^p(t) accelerates during the 
test, as the driver travelri to the catcher. achieves a constant 

value at the transient end point, when the driver hits the catcher. 
2.2.2 Aiialog Measurement Systems and Control Units 

Blade angles were measured Dy strain gages on the flapping 
flexures CFlgures 1 and 2). Similar flapping angle measurements are 
described extensively in References 27 and 29. The four strain gages 
on each flexure were wired into a Wheatstone bridge. The bridge 
remains balanced for blade edgewise and torsional deflections. The 
bridge responds only to blade flap^bending deflections. Each bridge 
circuit is separate. 3600 Hz 3 volt rms power is supplied, in the 
space fixed reference system, to each bridge, in the rotating 
reference system, through 3 of the circuit” in the 20 circuit slip 
ring assembly. This signal is amplified by an A. C. amplifier 
(Figure 9) and is recorded on a 6 channel F. M. analog tape recorder 
(Figure 9). The output level of the amplifiers was adjusted to give a 
one volt output signal per degree of blade flapping, when possible, so 
that flapping angle amplitudes could be easily monitored on oscillo- 
scopes during testing. The input to the flapping angle tape recording 
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channels could be switched for calibration purposes from the amplifiers 
to 2 D* C* voltage levels generated by a mercury battery circuit. 

There are control circuits for two variable speed motor drives 
(Figure 9). These controls are described in References 27 and 29. 

The control circuits of the inner shaft and rotor motor are wired so 
that a 2.8 volt D. C. signal is recorded on the tape recorder channels 
used for monitoring the shaft azimuth angles \mtil the motors are 
engaged. When the switch that engages a motor is thrown, the 2.8 volt 
signal Is asd.tched off and is replaced by the shaft azimuth angle 
measarement which varies between .2 and 1.8 volts. Thus the starting 
point of a transient test is clearly defined by a sudden voltage 
shift. For transient tests set 3, Table 4, the inner shaft motor 
switch also engages a solenoid that releases the inner shaft when the 
motor is engaged. This same switch was used in the spring excitation 
system to mark the start of the transient and release the driver arm 
(Figure 8 ). 

Rotor shaft position, 41 , and inner shaft position, were 
determined by blips generated by magnetic pick-ups (Figure 6 ) passing 
set screws in the shafts in harmonic testing. The signals generated 
by the pick-ups were recorded on two channels of the tape recorder 
(Figure 9). Rotor shaft and inner shaft travel at constant speed in 
haiiflonic testing, consequently one does not need continuous measurements 
of their position during a test. It is sufficient to know their rate 
and a position at one given time during a test. Their position at any 
other time in the test can be calculated from these two factors. 

During a transient test continuous measurements of the azimuttv^l 
position of the inner shaft, must be made so that the cyclic pitch 
input condition can be determined (see The blips cannot fulfill 

this function so tney were replaced. 

Two resolvers CFigure 7) driven through timing belt drives were 
installed to make continuous measurements of ^ and The resolvers 

were purcliased from a used parts catalog for a nominal sum. Operating 
specifications were not available. Experiments shoved that these 
resolvers produced a repeatable 2 per revolution signal instead of a 
repeatable 1 per revolution signal. The rotor resolver drive ratio was 
changed from 1:1 to 1:4* The inner shaft resolver drive ratio 
remained 1?1 since the starting position of the i ner shaft was known 
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at the beginning of the test and the repeatable 2 par revolution signal 
provided more accuracy In the measurement than a repeatable 1 per 
revolution signal could* Power was provided to the resolvers by audio 
oscillators operating at about 10,000 Hz* The output of the resolver 
is a 10,000 Hz signal attenuated by factors which vary with shaft 
position* The maximum modulation frequency Is 100 Hz* The resolver 
output passes through a full wave rectifier and then through a low 
pass filter with a frequency cut-off of approximately 1000 Hz* The 
resulting signal iB a repeatable 1 or 2 per shaft revolution waveform* 
la transient tests sets 1 and 2 (Table 4) the resolver signals were 
contaminated by inner resolver feedback through the single driving audio 
oscillator* This happened because the output Impedence of the 
oscillator was not low enough compared tc the lsq>edence of the two 
resolver circuits* Consequently there Is a 1 per revolution error sig- 
nal in the measurement of the Inner shaft azimuth position (4^^)* The 
error signal maximum amplitude is 5 percent of the maximum amplitude of 
the inner shaft resolver output. This error signal was eliminated in 
transient test set 3 by driving the two resolvers with different audio 
oscillators. 

The resolver signals were recorded on the analog tape recorder* 
Converting the resolver signals into ^ and measurements was done 
with digital algorithms In the conq>uter and is described in Section 3* 

Wake measurements were made with conventional hot wire anemometry. 
Two probes were calibrated in a calibration chamber for 0 to 20 ft/sec 
flow velocities* Linearization coefficients for the 4th order poly- 
nomial signal liaearizers were obtained from the calibrations* A 
probe stand (Figure 10), which allows probes to be positioned at 
various azimuthal and radial positions below the rotor, was constructed. 
Probe height can also be adjusted* The measurement stations used in 
the set of harmonic tests of Table 4 are shown in Figure 11* Wake 
measurement signals were recorded on the analog tape recorder* The wake 
measurement signals were also processed on a RMS voltmeter with a time 
period of 1 sec. and on a D. C. Voltmeter while the tests were in 
progress. 

The 6 channel Ampex SP-300 F* M. tape recorder (Figure 9) records 
D. C* voltage levels between -4 and -K volts* At the recording speed 
used for all tests, 3*75 ips, the recorded signal remains undistorted 
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up to frequeacies of 625 Hz. The recording level can be adjusted so 
that the signal actually recorded can be from .5 to 3»5 times the 
actual input signal. 1 to 1 input signal to recorded signal ratios 
requires very careful adjustment of the 6 level controls. The D. C. 
output level can be adjusted across a .5 volt range. At 3.75 Ips, the 
tape drive reaches steady state speed in about 3 Inches of tape from 
the stopped position. During a test, approximately 1 foot of tape is 
allowed to pass the record heads before recording is initiated. 

Measurements recorded on the tape recorder during harmonic tests 

were: 

channel 1 Rotor shaft blip for determining rotor azimuth position 

channel 2 Inner shaft blip for determining cyclic pitch azimuth 

position 
channel 3 ^ 

channel 4 Hot wire measurement of rotor wake. 

Measurements recorded on the tape recorder during transient tests 

were: 

channel 1 Rotor shaft resolver for continuous rotor azimuth position 

channel 2 Inner shaft resolver for continuously determining cyclic 

pitch azimutl. position 

channel 3 3^^ 

channel 4 3^ 

channel 5 3^ 

channel 6 3, 

4 

2.3 CALIBRATION AND TEST PROCEDURES 

In all tests the measurement signals were recorded on the 6 
channel F. M. anal<^g tape recorder. Some measurea^nt functions 
(Section 3) and all analysis (Section 5) were performed on minicomputers. 
The interface between analog signals and digital signals was provided 
by the analog-digital converters controlled by the computer. The 
testing procedure was designed around this interface. The analog to 
digital converter accepts signals betweer 0 and 2 volts. It assigns 
voltages of 0 to the integer value -512. It assigns voltages of 2 
to the integer value 512. All signal levels recorded on analog tape 
were maintained between 0 and 2 volts except for the pretransient 
signal recorded on the inner shaft azimuth position channel. The A/D 
converter can address 8 input lines and sample at a rate of once every 
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16 micro-seconds. 

2.3.1 Harmonic Test Calibration and Testing Procedure 
Calibration of the harmonic test flapping measurement required 

accurate knowledge of three calibration numbers: 

1. The aiq>lifler output voltage per inch of tip deflection^ which 
was measured by loading the blade tip and producing a known root 
moment and therefore a known tip deflection. 

The tape recorder recording level anqplification or attenuation 
which was carefully adjusted by recording known input signal 
levels and playing them back. 

3. The precise A/B converter mapping of voltage level into the 
±512 integer range, which was measured by sampling known D. C. 
voltage levels and displaying the results. The linearity of 
the conversion process was verified during this measurement. 

The hot wire wake measurements \«ere normalized so that 0 to 20 ft /sec 
yielded 0 to 2 volts output. The tape recorder record levels were 
carefully set and the A/D converter channels were studied for attenuation 
or amplification. The maximum signal level of the rotor and inner 
shaft blips have to be .3 to .8 volts for the program to find their 
relative position. 

The test procedure is as follows. The collective pitch and the 
ground height are set at the desired test condition. A gear ratio is 
chosen, for example one that yields ^ ^he hot wire probe is 

positioned at a particular radial position (Figure 11). The rotor 
is run up to speed. The tape recorder drive is engaged and a foot of 
tape allowed to pass the recording head then the recorder is put into 
record mode. Three feet of tape are recorded, a meter keeps track of 
tape position. While the signal is being recorded, the average wake 
value (the D. C. signal level) and the RMS wake value are recorded 
by hand. The rotor is stopped and the hot wire probe is moved to a new 
position. The sequence is repeated until 7 radial positions (Figure 11) 
of the hot wire probe have been recorded. The reversing gear is 
engaged so that = -.2 and a set of 7 wake measurements is made. 

This process is repeated until all 6 sets of gear ratios have been used. 
The whole test condition takes about 500 minutes (Table 3) to survey. 

2.3.2 Transient Test Calibration and Test Procedure 

Table 2 presents the interface between Instrumentation, mini- 
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computer, calibration procedure and testing procedure for the 
transient tests. Much greater use was made of the mini-computer in 
transient test calibration than in the frequency response test 
calibration. 

A set of calibration runs was conducted before each set of tests. 
Each blade was loaded at Its tip with a known mass which produced a 
known root moment and consequently a known flapping angle. The output 
of the amplifier for each blade flapping angle strain gage circuit was 
adjusted to give 1 volt per degree. The combined tape recorder and 
analog digital converter amplification and/or attenuation per channel 
were determined by recording 2 knoim D. C. voltage levels on (Table 2) 

4 flapping angle recorder channels. The tape was played back through 
the computer which averaged the large nusd>er of sampled results for the 
two voltage levels, and determined the ratio of tape recorder Input 
signal level to computer output signal level. The same digital volt- 
meter used to set the flapping angle amplifier output levels was used 
to measure the 2 tape recorder 0* C. voltage Input levels, consequently 
possible voltmeter calibration errors are avoided. There are then two 
sources of error In this calibration procedure; misplacement of the 
mass at the blade tips and D. C. drift of the F. M. tape recorder. 

The mass can be placed with 1% accuracy and the D. C. drift of the tape 
recorder yields a further 1% error source. Thus the calibration 
numbers are known to within 2% of their true values* 

The test procedure Is as follows. The collective pitch and the 
ground height are set at the desired test condition. The tape recorder 
Is engaged, the rotor motor is engaged - step D (Table 2), the Inner 
shaft motor is engaged - step E (Table 2), the tape Is stopped, and 
finally the rotor and Inner shaft motor are disengaged. The test takes 
less than 1 minute to perform. 

2.4 A COMPARISON OF TEST PROCEDURES AND TEST SET PHYSICAL CONDITIONS 

Table 3 presents a comparison of the time involved In conducting 
a test for the two testing methods. A test condition is defined as a 
combination of a particular value (first rotating flap-bending 
frequency), a particular collective pitch setting and a particular 
ground height. It takes about 1 minute to change the rotational speed 
and consequently the u)^ value, about 3 minutes to make a ground height 
change, and about 30 minutes to make a collective pitch change for both 
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Table 3 

Time Comparisons of Test Methods 


Harmonic Transient 

Test Test 


Time Required 

(til value 

1 min. 

1 min. 

to Cfaanga 
Test Condition 

Ground Height 

5 min. 

5 min. 

Collective Pitch 

30 min. 

30 min. 

Test Time 


500 min. 

1 min. 

Model Running 
Time 


30 min. 

1 min. 

Equivalent 
IBM 360 
Computer Time 
Required for 
Digitization 
and Analysis 

Harmonic Test 
Digitization and 
Fourier Series 

10 min* 


Transient Digitization 
and Fourier 
Coefficients of 1 D 


13 sec. 


Parameter Identification 


1-2 min. 

test procedures. 

The test time for a harmonic 

test is 500 

times that 


of a transient test. However, the model running time for the harmonic 
test is 80 times that of a transient test. Much time Is spent in the 
harmonic test changing gear ratios by hand. This task could be done 
automatically with a mechanically more sophisticated model, so the two 
procedures should be compared on the basis of running time rather than 
total testing time. 

Table 4 presents the physical condition of the rotor for each test 
set. For the early tests, the spread in collective pitch setting la 
large, however the average value of collective pitch for the rotor is 
quite accurate. In the case of the harmonic tests, the wake measure- 
ments take forms which are equivalent to averages over long periods, 
consequently individual blade differences are not Important. 

The measurement error of local rotor wake velocity does not 
exceed 4% of the true value. The error comes from Inevitable probe 
calibration error and local anemometer linearization error. The 
calibration error for blade flapping angle measurements does not 
exceed 4% of its true value for the harmonic test and does not 
exceed 1.5% for the transient test. This random error does not 
vary from measurement to measurement but is a constant error present 
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Table 4 


Physical Conditions of Tests 



average 

6 

0 

spread 
o{ 9^ 


ground 

height 

spread of 

pre-lag 

error 

spread of 

pre-cone 

error 

pitch 

error 

signal 

Harmonic 

Tests 

2° 

5° 

8 

i.7° 

1.18 

1.28 R 

.7° 

.3° 

5 X 

Transient 
Tests 
Set 1 

7 , 5 ^ 

t.4° 

1.18 

1.28 R 
1.02 R 
.78 R 

.7° 

.15° 

5 X 

Transient 
Tests 
Set 2 

3.3° 

6.3° 

9.3 

t.2° 

1.18 

1.28 R 
1.02 R 
.78 R 

.0° 

.0° 

5 X 

Transient 

5.0° 

i.i° 

1.18 

1.02 R 

.0° 

.0° 

5 X 

Tests 
Set 3 

5.0° 

i.i° 

1.25 

1.18 

1.12 

1.28 R 
1.02 R 
.78 R 

.0° 

.0° 

1 % 


in all measurements. The measurement error of the cyclic pitching 
excitation angle remains within 1.5% of its true value. The cyclic 
pitch input azimuth angle is known to within 5% for most tests and to 
within 1% for the last set of tests conducted. 
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3. DATA PROCESSING 


This chapter describes the manipulations of the transient data 
from the analog-digital Interface through the creation of digital data 
files In a form suitable for treatment by the Identification algorithm. 
Data processing of the harmonic test data Is adequately described In 
References 27, 28, and 29. 

3.1 THE MINI-CONPUTER SYSTEMS 

Data processing and the Implementation of the Identification 
algorithm were done on mlnl-coiq>uters . The work presented herein 
would not have been possible without the use of the modem mini- 
computer. The Identification algorithm Is essentially a digital 
algorithm which would cost too much to Implement In analog form. 

This work Is also not suitable for Implementation on a large scale 
conqputer for economic reasons. The conq>utatlonal cost of doing 
this work on a large scale system would have been approximately 
50 times what It cost to do the work on suitable mlnl-conqsuters . 

Tasks that require relatively small amounts of memoiry, large numbers 
of computations and few peripherals, are perfect for mini-computers 
and wasteful for large-scale systems since the task charge will 
include the overhead on large amounts of memory, a much more expensive 
processor, and many peripherals that are not necessary for the 
implementation of the task. The minis also offer the user the 
chance to Interact with his data much more frequently than a large- 
scale system would, thus offering great flexibility In implementation 
of data manipulation tasks and greater probability of early error 
detection. 

Program development was carried out on a PDP-8 which has a 24 K 
memory, two disk drives, and a fast line printer. The PDP-8 can 
communicate with a PDP-12. The POP-12 was used for data sampling 
of the analog tape and data reduction. The PDP-12 has an 8 K 
memory, two disk drives, two file oriented tape drives, analog- 
digital converters, a processor controlled scope with a 512x512 
matrix, and a teletype. The reduced data were written on tape and 
taken to another PDP-12, which has a hardware floating point processor 
and a 32 K memory, to be analyzed by the Identification algorithm. 

The results were then written on tape, taken back to the first 
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PDF- 12, sent to the PDP-8 and printed on the line-printer. 

3.2 THE DATA REDUCTION ALGORITHMS 

There are three groups of data reduction programs. They are: 

1.) sampling programs, 2.) calibration programs, and 3.) reduction 
programs. The sampling programs Interface with the tape operations 
of Table 2 as shown In Table 5. Calibration and reduction programs 
are applied to the resulting data sets to produce normalised data 
sets suitable for the parameter Identification algorithms (Table 5} . 
3.2.1 Sampling Programs 

Three F. M. tape sampling programs were written. These are 
assembly language programs that are similar to each other In many 
respects. They are written in assendily language because the Fortran 
callable subroutines that Initiate analog-digital conversion will 
not sample fast enough for these tests. 

The PDP-12 analog digital converter Is multiplexed to 16 Input 
channels. Eight of the input channel voltage levels are controlled 
by potentiometers set by the user. The remaining 8 are available for 
analog data Input. The PDP-12 control panel has 5 two-position sense 
switches which are controlled by the user. The states of these sense 
switches are program detectable and can therefore be used for 
branching purposes by the program. 

The three sampling programs are called RESOL, ARBSM, and SRSLA. 
They form the basic link between the tape recorder data and preliminary 
digital data files. The three programs have two basic operating modes: 
1.) sample and display, and 2.) sample and store. 

The sample and display modes of the three programs are Identical 
to one another. The sampling rate is controlled by the voltage level 
of one of the potentiometers. Two of the 6 data channels are sampled 
at a time. The digital values are scaled and then displayed on the 
PDP-12 scope. Two "zeroing lines which represent the middle of the 
channel input levels are displayed also. The D. C. output levels 
of the tape recorder can be adjusted so that the data signals are 
properly "zeroed" about the zeroing lines. Sense switches are used 
to choose between sets of 2 data channel displays. The sample and 
display mode allows the operator to make all necessary adjustment to 
the tape recorder output signals such as attenuation and D. C. level 
control so the the sanq>le and store modes will collect the proper data. 
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The program RESOL Is used exclusively for sampling resolver data 
for calibration purposes. RESOL is used on the taped data that lies 
between model manipulations A and B of Table 2. The user observes 
the output of tape channels 1 and 2 with RESOL's sample and display 
mode and moves directly to Che sample and score mode by changing the 
state of a sense switch when he determines that the correct Information 
is being displayed. The minimum possible time between a sample and 
store in memory on one channel and a sample and store on another 
channel is 20 microseconds. If the rotor velocity is 25 Hz , the 
rotor resolver signal Is periodic with a period of 40,000 microseconds 
and the pitch resolver is periodic with 20,000 microseconds. Taking 
360 samples/cycle of the rotor resolver (and therefore 180 samples/ 
cycle of the pitch resolver) yields 1 sample set every 111 microseconds. 
The rotor resolver sample and its respective pitch resolver sample 
are .2^apart because of the time b>>twcen 8anq>les. RESOL samples and 
stores the rotor resolver signal, then the pitch resolver signal, goes 
into an adjustable delay loop and repeats the process until 2047 samples 
of each signal are collected. After the sample and store is completed, 
the user can examine the sampled signal on the scope, by manipulating 
sense switches, to check signal Integrity. If the data is satisfactory, 
the user directs RESOL to write the sanple set onto an empty area of 
the disk in an integer formatted form. 

The program ARBSM is used for calibration data and pretransient 
data (Table 5). The user can observe tape channels 1 and 2, 3 and 4, 
or 5 and 6 through sense switch manipulations. When all input channels 
are adjusted to the user's satisfaction, and when the proper portion 
of tape is being read, the user triggers the sample and store mode 
by pressing a sense switch. ARBSM sequentially samples all 6 data 
channels, stores the results, goes into a delay loop and repeats 
until 1024 samples of each data channel are collected. The sampling 
rate is set at 60 sets of samples/cycle for a rotor veloctiy of 25 
cycles/sec. The angular difference between the first and last sample 
of each set of samples is less than 1°. In the analysis programs, it 
is assumed that each set of samples is collected simultaneously. This 
introduces a .3% error into the algorithm which is negligible. After 
the data is collected, the user can display it on the scope to check 
its validity (i.e. look for missing data channels). If it is 
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acceptable it la vrltten onto an empty area of the disk. 

The program SRSLA works exactly as ARBSH except that when the 
user exits "sanple and display” the program begins monitoring the 
pitch resolver data channel. When the signal level drops below 2 
volts on this channel, "sample and 8tore"ls initiated. The abrupt 
signal level drop signifies the beginning of the transient 
excitation. 

3.2.2 Calibration Programs 

The calibration program BTESTV is used after ARBSM (Table 5). 

It reads Che data left on Che disk by ARBSM and confutes Che average 
value of the sample sets for each data channel. It also gives the 
standard deviation of these signals. BTESTV is used only on D. C. 
calibration signals that are used to define flapping measurement 
attenuations. The standard deviations are always quite small - on 
Che order of .005 Volts - as compared Co Che difference between Che 
two calibration signals recorded on every channel - on the order of 
1. Volt. 

The calibration progreun CRSLT is used after RESOL Co create 
normalized pitch and rotor resolver Cables (Tables 2 and 5). CRSLT 
reads Che rotor resolver digitized data left on the disk by RESOL. 
The digitized waveform is approximately sinusoidal and is relatively 
noise free. As many of the 2047 available data points as possible 
should be used to create the most accurate Cable of signal aoqillCude 
versus azimuth angle. The program 'draws' an arbitrary line through 
the signal (Figure 12) of magnitude 0 (since the signal varies 
between -512 and 512). The signal is periodic with the number of 
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sample positions between the first crossing and the last odd numbered 
crossing divided by 1/2 the number of crossings - !• The signal is 
properly zeroed by averaging the maxima and minima (which are found 
by polynomial interpolation) between crossings and subtracting that 
nuinber from every point in the data net. The new maxima and minima 
are used to locally normalize the signal between crossings so that the 
signal is normalized between ±1. A 180 point/360 degree table of 
resolver magnitude versus azimuth angle is produced by averaging the 
normalized signal magnitudes for each azimuth angle. Linear inter- 
polation is used between data points. The resolver position is 
adjusted so that the zero azimuth angle occurs when the signal level 
is 0 with positive slope. 

The program CRSLT treats the pitch resolver data exactly the 
same way as the rotor resolver data. The end products of CRSLT are^ 
therefore, two 180 point /360 degree tables for both the rotor and 
pitch resolvers. 

It is very important that the pitch resolver signal level versus 
pitch azimuth angle be independent of pitch azimuth angle velocity. 

This was verified, in the frequency range of the tests, by recording 
pitch resolver signals of different frequencies, creating the 
corresponding normalized pitch resolver tables for each frequency, and 
then comparing the tables. Table values of like azimuth angle were 
within 1% of each other over a frequency range of .In to 1.5n. 

3.2.3 Reduction Programs and Diagnostics 

Each data set that is used for identification purposes is con- 
structed by combining a pre-transient and the corresponding transient 
data sets for a test condition. 

ARBSM is used to collect the pretransient data (Tables 2 and 5) . 
AVBETl produces normalized 60 point tables of each flapping signal 
magnitude versus azimuth angle as well as the average value of each 
signal frou the pretransient data. The ARBSM sampling rate is approxi- 
mately 60 points/rotor revolution or 1 sample every 6®. There are 
1024 samples or about 17 rotor revolutions of each data channel. Any 
one set of 60 points will not produce very accurate measurements 
of rotor azimuth angle or blade flapping magnitude. But since the 
waveforms are all periodic, the same sorts of data averaging 
techniques used to generate resolver calibration tables are used by 
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AVBETl to generate pretransient single blade data tables. The 
correct azimuth angles and the periodicity are determined from the 
rotor resolver data. The 4 single blade flapping measurement data 
sets are first zeroed to eliminate D. C. signal bias * using an 
average of each signal over the longest available Interval In which 
the signal Is periodic (these average values are also saved). The data 
sets are then averaged at each azimuth angle using linear Interpolation 
between data points. For Instance » at an azimuth angle of 36^ the table 
value consists of the average of 17 numbers. These numbers are the 
result of linear interpolation between the two nearest data points to 
the 36^azlmuth angle in each rotor revolution. Finally each table is 
normalized with the calibration numbers obtained with BTESTV. 

AVBETl takes the 4 single blade flapping tables and transforms 
them into 2 multiblade or forward and left tilt tables. Both the 
4 single blade and 2 multiblade tables are then Fourier analysed at 
frequencies of 1, 2» 3, and 4 times the rotor velocity. The Fourier 
coefficients were used as a diagnostic to make sure all blades were 
operating under the same condition. For exdiiq>le if the absolute 
value of the first Fourier coefficient of one blade was very much 
different from the other 3 blades, then that blade probably bad an 
incorrect collective pitch setting. If the differences were greater 
than 10% the test condition was thrown out. The absolute values 
the second Fourier coefficients were useful in diagnosing pre*lag 
and pre-cone setting errors (see Section 2*1.2). 

SRSLA is used to collect the transient data (Tables 2 and 5) 
that appears on the tape directly after the pretransient data. NORIDA 
produces normalized pitch resolver azimuth angle measurements and 
the 4 single blade flapping angle measurements. Of the 6 data tracks 
only the rotor resolver signal is periodic. This signal is used to 
determine the number of data points per rotor revolution and the posi- 
tion of the first zero azimuth angle (the azimuth angle of any data 
point can be calculated using these two quantities) using the averaging 
techniques that were used to create the resolver calibration tables. 

The pitch resolver signal is zeroed and then normalized about local 
maxima and minima (Figure 13). The data values are then compareo to 
the pitch resolver normalized amplitude versus pitch azimuth angle 
calibration table point by point and are converted into pitch azimuth 
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local minima 


-512 

0 sample position 400 

Figure 12 Transient Pitcli Reaolver Sample Window 


angle measurements. The 4 single blade measurement sets are left 
unchanged by NORIDA. 

HVOHUL uses the normalized pri'translent data set and the data 
set created by NORIDA (Table 5) to create the data set used for 
identification purposes. The identification data set has 5 
measurements: 

1. ) Number of points per cycle 

2. ) Point of first zero rotor azimuth angle 

3. ) Pitch azimuth angle 

4. ) Left tilt 

5. ) Forward tilt. 

The data set has these properties: 

1. ) Measurements 3) > 4) > and 5) start with approximately zero 

value 

2. ) The tilts are normalized with the amplitude of the excitation 

eccentric. 

The rationale for using these measurements with these properties is 
given in Chapter 5. 

HVOMUL uses the averages of the pretransient single blade measure- 
ments to eliminate transient single blade D. C. signal measurement 
bias. This is possible since the actual signals are recorded only 
seconds apart (Table 2) which allows little time for instrumentation 
D. C. drift. The 4 single blade transient data sets are converted 
to the 2 multiblade measurements using the azimuth angle measurement 
from NORIDA. These multiblade measurements are 'zeroed' with the 
pretransient values of forward and left tilt (which are also the 
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Initial conditions of the forward and left tilt for the transient 
tests) . The pitch azimuth angle is . ' zeroed* with the average of the 
first 10 samples of Itself. This Is possible because the transient 
starts about 40 samples Into the data set. 
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4. AHALYTICAL TREATMENT OF THE HOVERING UNSTEADY ROTOR WAKE 

COUPLING 


The model rotor Is essentially infinitely stiff in the chordwise 
and torsional directions for the frequency range of pitching excitation 
that the model is subjected to (Table 2). Therefore only rotor flap- 
wake interactions need be analysed to predict the rotor model response. 
For a N bladed rotor operating at zero advance ratio, assuming 
straight blades flexibly binged at the rotor center, the single 
blade equations of motion in a rotating reference system are. Refer- 
ence 11, 


. ) * • * y I 


(4) 


By 

vliere A»-^S A** Is a function of the wake distribution over the radius, 
and is the generalized wake coordinate. The actual wake velocity is 
obtained by multiplying with the appropriate radial wake distri- 
bution. 

At this point, no assumptions have been made as to the form or 
structure cf the wake. The quantities A*’Xj^ represent the sum of 
the effects of the wake on the flapping motions of the rotor. In 
Reference 11 it is established that the hub moments produced by the 
flapping restraining forces (inertia and spring forces) about the 
pitching and rolling axes through the rotor center are proportional 
to a first order expression of the generalized induced flow coordinate 


+ TVj = -L{^ cosipj^ } 

k“ L 

b .. 2 

''ii ^ 

k=l 


(5) 


where L is the induced flow gain factor and x the time constant 
describing the delayed build-up of the wake, and are the 
generalized multiblade induced flow coordinates and are hereafter 
referred to as the dynamic wake. Transforming the single blade 
variables to multiblade variables using the transformations 
(Reference 35) 
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®k “ " ■*“ ®ii‘^°®'hc 

\ ■ '’i“® *k * "11°*°^ 


yields 

3j. + A3j + (u>J - l)3j. + 23jj. + A6jj - A(6jj + v^) a) 

-23j. - A3j. + ‘gj-j + ASjj. + (J[ -l)3jj - A(-0j + Vj.j) b) 

2 . (7 

Vj + TVj = -L {Bj + (oil -DBj + 2Bjj.} c) 

Vii + TVJ.J. =■ -L (Bjj + (a.2 -DBjj + 2 Bj} d) 

if one assumes that the d}niaiiiic wake Is the only part of 

the wake that is important to flapping dynamics. Rearranging (7) 
yields multiblade equations in state variable form x»Cx+Bu, 
that is 



u 130 / ^ l^ T* AL . _ 1+AL 
where P2=(w^-1), L*« — , x* = — — . 

No assumptions have been made as to the numerical values of L 

and T, which define the dynamic wake. Only a mathematical framework 

for their use has been provided. Given this framework and a proper 

set of experiments, the parameter identification algorithm developed 

in Section 5 can be used to find numerical values of L and x for the 

rotor upon which the experiments were conducted. This mathematical 

framework has one very interesting feature. If the rotor is excited 

2 

at a constant frequency to « in pitch, then as can be seen 

in (5), there can be no dynamic wake (v^ and are 0). This 
feature of the mathematical model can be easily checked, at least 
qualitatively, by performing a simple experiment with the model rotor. 
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A hot wire probe was placed at the .02R blade station .2R 
below the rotor. The rotor was excited with a pitch stirring fre- 
quency ofu al.2or(i>B«2. The rotor speed was varied from low to 
high rpm so that varied from 1.35 to 1.1. At a 1.20, the 
dynamic wake should disappear. Figure 12 presents a series of 
oscilloscope photographs of the hot wire output versus w^. In the 
upper trace, each vertical mark represents one rotor revolution. 

The lower trace is the hot wire measurement of the wake. The 
"dynamic" part of the wake, if it exists, should be periodic with 
5 rotor revolutions. The photographs show that there is indeed a 
wake component periodic with 5 rotor revolutions that goes toward 0 
as goes to 1.2. Thus the mathematical framework is promising in 
that it predicts a measured wake phenomena. 

Some Interesting variations on and modifications of this 
mathematical framework are now considered. Rotor drag should 
induce a flow in the direction of rotor rotation, therefore creating 


feedback between and Thus the right hand sides of equations 

(7), c) and d) become 

TV + V + Uy 

I + I 1- HVjj 


+ ''ll 


Elements 5»6 and 6,5 of the C matrix in equation (8) become 

u 

— and -H’ respectively. 

A logical extension of the time delayed Induced flow Is to add 


an acceleration term. The right hand sides of (?), c) and d) become 
TiVt + T 2 VX 

>*• ^ ^ /srkV 


+ '‘’2V11 + Vii 


In state variable form this model becomes 
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where T 



1+LA 

> 



and L' 



If one assumes that the wake gain factor L la adequate for 
modeling purposes, that t « 0, equation (8) reduces to 



( 12 ) 


A 

where A* =» similar to the Curtis and Shape y » 22. 

This Is a quasi-steady wake model where the wake Instantly responds 
to rotor motions. 

Equations (12), (8), and (11), In order of Increasing complexity, 
represent three possible mathematical frameworks for treating 
dynamic rotor wake coupling. They are, respectively, quaslateady, 
first order unsteady, and second order unsteady treatments of the 
wake. Equation (9) introduced an interesting variation on (8). The 
applicability of each of these constructions can be easily studied 
and compared using the experiments described In Section 2 and the 
parameter identification algorithm developed in Section 5. 

A simple momentum theory o* dynamic rotor Induced flow is presented 
in Reference 13. It is based on the assumption that the axial Induced 
flow is uniform over the disk and that the harmonic Induced flow given 
in Equation (5) varies linearly with the radial distance from the rotor 
center. The air mass participating In the angular pitching and rolling 
motion is determined from potential flow theory for a solid circular 
disk. Under these assumptions and under the further assumption of small 
linear perturbations about trim one obtains for zero advance ratio the 
following expressions for L and x: 


L ■ aa/C2y v^) , x « .113/ (13) 

Instead of L, Reference 18 gives the factor of the rotor moment coef- 
ficients as while Equation (8) is written in terms of blade 

moments divided by the blade moment of inertia, To obtain the 

rotor moment coefficient one must multiply M^^/I by I(b/2)/pnR^, With 
I » apcR^/y and o » bc/TiR this factor becomes aa/2y and the expression 
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for L In Equation (13) follows. With a ■ 5.6, o a .154, -y ■ 3.8 we 
have ao/2Y “ .113. Thus, according to the theory 

L - T - .113/ (14) 

Theoretical values of v can be calculated from blade element 

o 

momentum theory. The experiments were conducted In the vicinity of 
the ground so the momentum theory must Include the effect of the 
ground on the Induced Inflow. There are two conflicting methods of 
doing this. A theory that Includes the ground height effect on 
Induced flow for constant circulation (twisted) rotors Is described In 
Reference 36. The theory conflicts with the results of a comprehensive 
set of tests which measured the hovering performance of various helicop- 
ters as a function of ground height. This set of tests Is described 
In Reference 37. 

A Is a function of the tip loss factor B (Reference 32). B 
decreases with Increasing collective pitch and/or induced Inflow. 

Therefore A can be expected to decrease with increasing collective 
pitch. 
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5. DEVELOPMENT OF THE IDENTIFICATION ALGORITHM 

The proposed perturbation models of Section 4 have from 6 to 
8 state variables and from 4 to 6 parameters. The Identification 
method chosen had to be able to Identify all parameters from the 
transient cyclic pitch stirring tests so that no theoretical assump- 
tions beyond the form of the equations need be made. Only 2 state 
variables were measured, the left and forward flap-bending deflections. 
Therefore only two of the 6 to 8 state variables are directly 
related to the measurements. Thus the Identification method must 
not only Identify the parameters, but must also estimate the 
remaining non-measured state variables, the flapping rates and Inflow 
states. In other words, the method must Identify the main system 
parameters P and A and the feedback parameters without measuring a 
feedback signal and without measuring half of the main system state 
variables. In control system engineering this would be considered a 
formidable task tliat Is solved here by application of the proper 
Identification algorithm to the problem along with some rather 
stringent requirements on the type and duration of the Input transient. 

There are many parameter Identification methods available In 
the literature. Extensive simulation work was done to pick out the 
best method for identifying the parameter models described In 
Section 4 before the experiments were performed. This simulation 
work Is described In Refei ces 38 and 39. The best method was 
found to be a modification of the complete maximum likelihood 
Identification algorithm described in Reference 40. 

The complete maximum likelihood Identification algorithm can be 
used to Identify system parameters In the presence of both process 
and measurement noise. In order to apply the complete maximum 
Identification algorithm a "likelihood function" has to be maximized 
subject to several non-linear constra. nt equations (Reference 40). 

This is a formidable task. The likelihood 'function can be 
maximized using a Newton-Raphsou Iterative technique if one assumes 
ttiat the best state estimate covariances are constant throughout 
the data set and that the Kalman filter gain is also constant 
throughout the data s'^t. If one further assumes that there Is no 
process noise, then the Kalman filter gain is zero and the filtering 
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equations reduce to the system equations. The following Identification 
algorithm Is based on the foregoing assumptions. 

5.1 THEORETICAL FORMULATION 


The Iterative method developed assumes sytcem and measurement 
equations In the form 

X « f(x>u,d) + w (15) 

y ■ h(x) + V (16) 

where x Is the state vector, u the input vector, 6 the parameter 
vector, w the process noise vector, y the measurement vector, and v 
the measurement noise vector. The parameter vector 8 Is estimated 
In such a way that the criteria function or cost function 


T -1 

J(0) « (1/2)53 '^4 ® minimum 

jol J J 


(17) 


where 

V. 


yj - h(*j) 


(18) 


Here Is called the Innovation vector at time J 
(see Reference 40), is the best estimate of the state at time J, N 
is the nundier of measurements, and B is the measurement equation error 
covariance matrix and was updated before proceeding from the k to 
the k+1 iteration by 


J. T 

\+i “ 23 (y-\(x))j (y-\(x))j. 


(19) 
.th 


B in this form is a direct measure of the modeling error of the k 
iteration and it is convenient to think of B as the "modeling error." 
In the complete method B is a function of J. 

The cost function (17) was minimized using what is called quasi- 
linearization with the modified Newton-Raphson method (see Reference 
40). This method makes use of the information or sensitivity matrix, 
M, where 


M ^ 32j/302 . (20) 

M ^ gives the Cramer-Rao lower bound for the parameter covariances 
and has a very special Importance that will be demonstrated later on. 
The problem is to find a zero of the gradient of the cost function 
3J/ 30 B 0. Taking the truncated two term Taylor series expansion of 
the gradient about the value of 0 in the k** iteration. 
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OJ/36)j^^^ - OJ/30)j^ + O2 j/302)^ - 0 (21) 

one obtains with Equation (.20; the parameter update Increment >o be 
used for the k-t*! Iteration 


-k+i ^'Vi • -V<w**)k 

In order to obtain the parameter update needs the first 

and second partlals of J(6), 

^ T -1 

?J/30 - £ (av/36); B ^ V, (23) 

j-1 ’ J 

N 


48 . 




( 22 ) 


M 


)2j/ 302 - £ (3v/90)^ b“^ (3v/30). 

j-1 ^ ^ 


(24) 


where higher order partlals of are neglected (quasl-llnearlzatlon) . 
Computationally, this means that Vj and 9v^/30 are required, v Is 
obtained from (18) where x Is obtained from solving the system 
equation 

s 

X - f(x,6,u) (25) 


and 0 Is the best estimate of the parameters from the previous 
iteration. Now 


3v/30 = -3h/3x| - 3x/30 (26) 

where 3x/d6 Is obtained from solving the sensitivity equations 

3x/30 = 3f/30 + 3f/3x|^^- 3x/30 (27) 

The Initial conditions for dx/db are zero except when x(0) is 
identified as part of the parameters 6* In this case the initial 
partlals have the value one, 

A closer look at (24) and (26) will show that M is not a function 
of the innovation sequence but is only a function of the theoretical 
model and the input u weighted with the "modeling error" B from the 
previous iteration. The diagonal terms of the inverted sensitivity 
matrix M ^ are theoretically a lower bound for estimates of the 
parameter errors. As the theoretical model approaches a true 
representation of the experimental model the diagonal teras become a 
better and better lower bound. 

So computationally what is required is: 
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1) Make reasonable Initial estimates of parameter values and 

2) Solve simultaneously system and sensitivity equations. 

3) CoDq>ute the Innovation sequence* 3 J/dd, H, and 

4) Invert M and calculate parameter Iteration. 

5) Use updated parameters in step 1). 

6) Bepeat the process until the parameters converge. 

5.2 PRACTICAL COMSIDEBATIONS 

The above theoretical formulation of the Identification algorithm 
must be properly implemented. Tools that assure the proper applica** 
tlon of (his algorithm to data sets must be developed. These tools 
must show the user that the algorithm Is doing i(.s Job of identifying 
parameters and that the Identification process Is working properly. 

Convergence of parameters with successive parameter Iterations 
Is the most obvious place to start when comparing paraneter identifi- 
cation methods and when determining which Inputs excite the model 
sufficiently to identify parameters. If parameter values fall to 
converge then either the method is at fault or the input Is 
Insufficient. 

The square root of the average of the "modeling error" or B matrix 
diagonal 

RR = (28) 

hereafter called the "fit factor" » is a convenient direct measure of 
how veil the theoretical model fits the experimental data set. 

Determining the "best" input function for identifying parameters 
is the most difficult task that faces the experimenter who wishes to 
use parameter idantificatlon methods. Several works on "optimal" 
input design have been written (References 40 and 41) • The 
experimenter usually does not have the luxury of being able to use 
a theoretically optimal input in his experiment because of practical 
considerations such as physical limitations on his model, But» the 
experimenter can usually devise several different kinds of inputs with 
which to excite his experimental model. Therefore, he needs methods 
of comparing tne results obtained with various inputs. The experimenter 
also has to choose how much of a data set to use in the identification 
process* If he does not use enough data his Identification algorithm 
might not converge. ' Or worse, the theoretical parameters might be 
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relatively insensitive to the input function up to that point so that 
even if the parameter iterations converge, the coat function is 
locally flat in the area of the Identified parameter and a wide 
range of parameter values would give approximately the same fit 
between experiment and theory. If be uses data beyond the point 
where the experiment is generating new information that is relevant 
to the theoretical model, then he is using excessive processing time 
and is taking the risk of biasing his parameter estimates with 
Irrelevant Information. 

Convergence and fit factors can be used to compare parameter 
Identification methods, theoretical models, and inputs, but these 
two tools by themselves are Inadequate for answering all the 
experimenter's questions with regard to his data and models. A 
third tool that the experimenter might use that can answer his 
questions involves conq>utation of the cost criteria, J, as a function 
of parameter value and data length. This Involves excessive conq^uta- 
tlon but can be used as a last resort. 

The proper use ot M ^ is the tool the experimenter needs to 
answer bis questions about the relative value of input functions and 
how much data from his data sets he needs to use. The use of M ^ for 
determining proper data length is described in References 38 and 42. 

M ^ is a monotonically decreasing function with time (or data length) 
and is a function of the theoretical model and the experimental 
input weighted by a constant (the modeling error matrix); see 
Equation ( 24) . Section 5.1. The experimental input "drives" the 
value of M ^ toward zero. Inputs can be compared by observing how 
fast they drive M ^ toward zero and ultimately by how near M ^ gets to 
zero. If an element of M ^ is flat in a particular section of the 
data, then the algorithm is not making use of the information in 
that section of data to lnq>rove Che estimate of the parameter asso- 
ciated with that element of M . If there is a particular data 
length where all the elements of M are flat and remain flat for the 
rest of the available data length, Chen no further information is 
gained by using data past that particular data length. This 
establishes a clearly defined limit on how much data one has Co 
analyse to identify a particular theoretical model. Frequently in 
early simulation and experimental work insufficient amounts of data 
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were used for parameter identification and the unsatisfactory results 
were misinterpreted. The time history of M provides a tool for 
avoiding this pitfall. 

Caution should be used when applying the M~ criteria. The 
''criteria" Is only as good as the theoretical model. It "approaches" 
reality only as the theoretical model "approaches" a true model of the 
experiment. It should only be used after the parameters have con- 
verged to a constant value and therefore offer the best model that 
the theoretical framework can supply. 

Convergence, fit factor, and M studies form a powerful set 
of tools for using the parameter identification algorithm. These 
three tools are applied intensively to the data sets that were ^ 

generated in the experiments described in Chapter 2. 

Verification of the identified model is acconq>llshed by seeing 
how well the model response predicts experimental responses not 
used for identification purposes. 

5.3 IMPLEMENTATION OF THE ALGORir 1 

Blindly implementing the theoretical algorithm as stated in 
Section S.l would lead to an unnecessarily long and cooq>licated 
computer program. Careful consideration during implementation can 
greatly simplify the resulting program and even iiq>rove the informa- 
tion content of the progr 2 un output to the user. 

The algorithm was Implemented using Fortran IV in 16 K of core 
memory of the 32 K PDP-12 described in Section 3.1. The program is 
in Appendix 9.2. The prime consideration in implementing the 
algorithm was to use as little core memory as possible so that the 
program could be rmi on as many machines as possible. The secondary 
consideration was to minimize the number of computations in the 
algorithm. 

5.3.1 Measurements 

The physical measurements taken are (Section 2.2.2) the blade 
flapping angles in the rotating coordinate system ^ 2 * @ 3 > 3^* 

The related state variables are the blade flapping angles in the fixed 
coordinate system, 3^- and 3^^* The true measurement equations are 
then (Equation (16), Section 5.1) 

3j^ » 3j costjij^ + 3jj sintpj^ for k=l,2,3,4 (29) 

m 
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Instead of using C29) as the measurement equation, the coordinate 
transformation 

3t ® ((^1 “ )*cosi(( - (3, - 3, )*sinij;) / 2 

^m •‘•m Jm *m ‘*m (30) 

6tt = CC6, - 30*®i“’l' + (6, - 3, )*cos(J») / 2 

*** 1)11 

which is the Inverse of (6), was performed during the preprocessing 
of the data sets (^action 3.2*3) so that the measurement equation 
becomes 






«II *2 


( 31 ) 


This ±8 done without loss of information. 

Using (31) as the measurement equation cuts the measurement 
storage requirements in half since there are only half as many multi*- 
blade measurements as single blade measurements. Each data set 
consists of approximately 400 measurement vectors. Each number 
consumes 3 words of storage. Therefore, using multiblade measurement 
vectors saves 1.6 K of storage. If (29) were used as the measurement 
equation then the best estimates of 3^ and 3^^ from the system 
equation would have to be translated by (29) into best estimates of 
the 3|^'s for each iteration. Using (31) as the measurement equation 
eliminates this transformation for each Iteration. Thus using (3l) 
as the measurement equation saves storage and saves computations. 
5.3.2 Choice of Parameters 

The parameters that have direct physical significance appear in 
2 

Equation (7) as (o)i-l). A, L, and x. When (7) is written in the 
state variable form, (8), the physically significant parameters 
become lumped into the groupings (o)^-l) * P2, A, L* * AL/x, and 
X* s (l+AL)/x. The parameter identification literature (Reference 
40) recommends that lumped parameters be used in the identification 
algorithm. All simulation work was performed identifying lumped 
parameters (References 38 and 39) • Initially the algorithm applied 
to real data sets was programmed to identify the lumped parameters. 
However it was found that by identifying the physically significant 
parameters the algorithm converged more rapidly than lumped para- 
meter identification for real data sets. For example see Table 6. 
Identifying physical parameters requires using about 10 percent more 
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Table 6 

Parameter Choice Method Studies 


Iteration A L t 3^(0) 3jj(0) Vj(0) Vjj(O) RR 

I Physical Parameters plus bias* steady pitch stirring trim 



0 

.450 

6.0 

8.0 

.0 

.0 




.079 


1 

.283 

2.0 

4.5 

-.05 

-.03 




.037 


2 

.303 

2.6 

6.4 

-.05 

-.02 




.034 


3 

.307 

2.6 

6.2 

-.05 

-.02 




.034 


4 

.307 

2.6 

6.2 

-.05 

-.02 




.034 


5 

.307 

2.6 

6.2 

-.05 

-.02 





II 

Lumped Parameters 

plus bias 







0 

.450 

6.0 

8.0 

.0 

.0 




.079 


1 

.281 

2.0 

5.6 

-.05 

-.02 




.038 


2 

.302 

2.7 

6.6 

-.05 

-.02 




.036 


3 

.307 

2.6 

6.2 

-.05 

-.02 




.034 


4 

.307 

2.6 

6.2 

-.05 

-.02 




.034 


5 

.307 

2.6 

6.2 

-.05 

-.02 





III 

Physical Parameters plus bias and initial conditions, 0 

trim 



0 

.450 

6.0 

8.0 

.0 

.0 -.19 

.33 

.0 

.0 

.092 


1 

.282 

1.2 

3.2 

.01 

.02 -.19 

.31 

.33 

-.23 

.039 


2 

.293 

2.3 

5.8 

.00 

.03 -.18 

.31 

.42 

-.33 

.028 


3 

.302 

2.4 

5.3 

,00 

.03 -.18 

.30 

.51 

-.34 

.027 


4 

.300 

2.4 

5.8 

.00 

.03 -.18 

.30 

.48 

-.35 

.027 


5 

.302 

2.4 

5.7 

.00 

.03 -.18 

.30 

.49 

-.35 



Test Condition: 0^ = 9.3°, h/R =1.28 

Computer Time (PDP-12 with FPP) 

per iteration iterations required total time 

for complete convergence required 


I 

1. 

minute 

3 

3. 

minutes 

II 

1. 

minute 

4 

4. 

minutes 

III 

3. 

minutes 

«7 

21. 

minutes 



so 


terms in the sensitivity equations than Identifying luiiq>ed parameters, 
but It can save 1 to 2 Iterations. Therefore If It takes S Iterations 
for the lusq>ed parameter method to converge. It will take only 3 or 4 
Iterations for the physical parameter method to converge, using 70 to 
90 percent of the computer time of the lunqped parameter method. The 
great advantage that Identifying physical parameters over luaqied 
parameters has Is that the matrix has direct physical significance 
with respect to the physical parameters. The elements of M ^ can be 
used to directly determine the relative Importance of the physical 
parameters to the theoretical model for a particular Input. 

5.3.3 Sensitivity Equations and Matrix 

The Identification algorithm requires the solution of one 
sensitivity equation, (27), of the same order as the state equation 
for each parao^ter identified. For Instance, using the model of 
equation (8) and Identifying the parameter L yields the sensitivity 
equation (32): 
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The 6x6 matrix of (32) is identical to the 6x6 matrix of Equation 
(8). The sensitivity equation C32) has the same form 
as the system equation (8), with a new forcing function, which Is 
3f/3x|^-3x/3L or the second part of (3^. All of the sensitivity 
equations are of this form, which makes them very easy to program. 
Using a DO loop, the system equation can be used once for the 


(32) 
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system equation» once for each parameter to be identified, and then 
the forcing function terms can be added outside the DO loop at each 
time step. This greatly reduces programming time. 

In order to perform the parameter iteration the sums of 
Equations (23) and (24) must be calculated and M must be inverted. 

If one calculates these sums and then does the Inversion at the end 
of each iteration, then one needs to store values equal to the length 
of the data set used times the dimension of the measurement vector 
times 1 plus the number of parameters identified. For example if one 
uses all 400 measurement vectors in the typical data set and identi- 
fies 4 parameters then 400 x 2 x (1+4) « 4000 values must be stored, 
which requires 12 K words of computer storage. Instead of calculating 
these sums at the end of each iteration, partial sums were calculated 
after each 15 time steps (or data vectors), thus reducing the 
storage requirement to 15 x 2 x (1+4) = 150 values or .45 K words. 

By inverting the sum of the partial sums which make up M after each 
15 time steps one can readily obtain the time history of M which 
one wants on certain occasions for diagnostic purposes (Section 5.2) • 
This capability is incorporated into the program and can be switched 
on or off at will. 

5.3.4 Initial Conditions 

Whan not all state variables are measurable, as is the c^se here, 
the identification of the initial state vector in addition to the 
system parameters is usually required. The model presented in 
Equation (8) has 6 initial states and 4 system parameters. Therefore 
if initial conditions must be identified, 10 sensitivity equations are 
required . 

The transient excitations are designed so that up to t*0, 6^ and 
6^^ are in a steady-state condition. Equations (2) and (3). There- 
fore the multiblade coordinates of the models must have reached a 

a 

steady state condition, so the Initial conditions of rates 3^.(0) amd 
3^^(0) of (8) must be zero. This eliminates 2 Initial condition 
parameter Identifications and leads to 4 sensitivity equations for 
Initial condition parameters. 

The fact that both Inputs (6^ and 0^^.) and non-rate states 
(3i> Vj, and In steady state at t~0 can be used to 

eliminate the rest of the Initial condition parameters. One simply 
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''perturbs" the model about the steady state condition. The new 

measurements are then the difference between their old measured 

value and their pre-trenslent, or equivalently their tsO, value 

(3, ® &T ~^T ^TT " ^TT “^TT CO)). This requires subtracting 

im ^iBj. iim lim^ ^^®pt 

a constant vector from every measurement vector in the data set. 

This procedure is described in Section 3.2.3. The new inputs 

become the difference between their old measured value and their 

pre- transient, or equivalently their t»0, value (0_ ° 6--0_(0), 

I It Ipt 

6-- « 8-t-6ttCO)). Whrn the measurements and inputs are written in 

11 iij. iipt 

this way, all inputs up to t»0 are zero, therefore all of the 
perturbation model states up to t»0 are also zero. Now all states 
have zero initial conditions and no initial condition parameters 
need be identified. Table 6 presents a conparlson of the parameters 
identified from a typical data set both with and without initial 
condition identification where the appropriate forms of the input 
and measurement were used in both cases. There are no essential 
differences between the results. 

Careful design of the input transients reduced the 
number of potential parameters in Equation (8) from 10 to 8. 

Careful selection of the form of the measurement vector and input 
vector further reduced the potential parameters from 8 to 4. 
Potential computer tine was more than halved and the complexity of 
the problem was greatly reduced. 

5.3.5 Measurement Bias 

A measurement bias is an arbitrary constant value present in 
every measurement of the data set. The bias can be considered to 
be a parameter and can be identified using the identification 
algorithm. There are at least two ways of including measurement 
bias in the algorithm. 

In the method found in most parameter identification litera-* 
ture, the bias is handled in the system equations. New state 
variables are defined for the measured states. For exanq>le one 
would make the substitutions 

h ‘ h 

old ^new 


( 33 ) 
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In Equation (8). The equation would remain the same except for 


the new 

terms 
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(34) 


which would have to be added to the right hand side. The bias 
and 7^^ appear explicitly in the system equation. The cost 
function remains the same. Writing that cost function using slightly 
different notation than (17) gives 


JC9) - % lie, - 05) 

where B is the measured state vector and B the best estimate 
m new 

of that measured state vector* The bias appears Implicitly In the 
cost function, therefore finding dJ/3'F^ and requires 

solving two more sensitivity equations* In general, one must solve 
a new sensitivity equation of the order of the state equation for 
each bias identified using this method* 

There is no reason not to consider bias explicitly in the cost 
function, which then becomes 


J(6) - % l|8„- 5„1 j-^IIb 06) 

where *3 is the bias vector* Since the bias only appears explicitly, 
aJ/a7j and 3J/a7j-j are Just ~(3j^“8j) and “(Sjj- ~3j.j.) which 
are already available. Thus no sensitivity equations need be 
solved In order to identify bias. Indeed the quantities one needs 
to conq>ute to identify bias have already been computed. 

The two methods ere mathematically identical « but the second 
is conqiutationally much more efficient. 

5.4 PBOGKAM DESCRIPTION 

A partial listing of the program used to identify the para- 
meters of the model described by Equation (8) is in Appendix 9.2. 

The program will not be described in great detail here because it 
is Internally documented. 

The differential equations were solved forward in time using 
a modification of a standard 4th order Runga-Kutta subroutine from 
the IBM Scientific Subroutine Package. The modification was the 
elimination of the provision for halving the time Interval of 
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Integration in tlu event of a computed order of error becoming larger 
than a user specified limit. Previous experience has shown that the 
program almost never made use of the interval halving procedure. 
Elimination of this procedure reduced the subroutine size by one 
third. 

The matrices were inverted using a standard matrix inversion 
subroutine OHNV) from the IBM Scientific Subroutine Package. 
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6. RESULTS AMD DISCUSSION 

The harmonic test results are presented in Section 6.1. Data 
length studies, idilch show how much data has to be used for proper 
parameter Identification, are presented In Section 6.2. Identified 
parameters are presented in Section 6.3. In Section 6.4 prediction 
studies that show how well, or poorly, the Identified models predict 
transient response tests are presented. There are also some prediction 
results presented In Section 6.1 which show how Identified models predict 
frequency response tests. Those prec-ictlons presented In Section 6.1 are 
discussed In Section 6.4. Figures are presented at the end of each 
section. 

6.1 HABMONIC TESTS 

The complete set of the harmonic test (Table 4) results, con- 
ducted as described in Section 2.3.1, are presented In Reference 34. 
Sanq>le results of the wake measurements made with a hat wire anemometer 
for the and (i>j^~1.17 case are presented In Figures 14a, b, c and d. 

Figures 14a and b present two views of the measurements of the 
time average of the rotor downwash velocity, U. In Figure 14a V Is 
plotted versus pitch stirring excitation frequency at different radial 
positions. This plot shows that the average wake Is largely Inde- 
pendent of the pitch stirring frequency of excitation. The Mf 

data (Ref erency 34) shows slightly less frequency dependence. W versus 
radial position at different pitch stirring excitation frequencies Is 
plotted In Figure 14b. The maximum value of ¥ always occurs at 
approximately r/R ° .8. Since the wake below a rotor rapidly contracts 
and since the measurements were taken .12R below the rotor this Is a 
reasonable result. 

The ¥ data at the three collective pitch settings were averaged 
over frequency at each radial position. A mean wake velocity at each 
6^ was determined by an r weighted average of the frequency averaged 
radial measurement values. These values are presented In Table 7 
where the results are given In units of radians. 
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Table 7 



Average Inflow at 

h/R * 1.28 


®o 

^0 measured 

'’036 

'’037 

L *• T 

2° 

• Glu 

.014 

.019 

7.0 

50 

.022 

.029 

.040 

5.1 

80 

.033 

.040 

.057 

3.4 


Theoretical values of Induced Inflow were calculated as a function of 


radius and 6^ using blade element momentum theory at h/R ■ The 
theoretical mean wake velocity at each 8^ was determined by an r 
weighted average of the theoretical radial wake velocities. Two 
methods of accounting for the effect of ground on the average Induced 
Inflow were used. Reference 32 used the results of Reference 36 to 
compute plots of h/R and blade radial position versus the reduction of 
Induced Inflow for a constant circulation rotor (Section 4). The 
reduction factor at r .7R was used to calculate the values shown 

In Table 7 from the theoretical values of average Inflow at h/R ■■ ». 
Reference 37 derived an empirical formula from tests for computing the 
reduction of average Inflow from h/R (Section 4). This reduction 
factor was used to calculate the values shown In Table 7. The 

L s T values In Table 7 were calculated with Equation 14 using the 

. values. Since the values agree better with the tests 

0 measured 036 ^ 

than the values, theoretical values of L and t from (14) were 
calculated using values and are presented In Section 6.3. 

Figures 14c and d present two views of the measurements of the 
absolute value of the first Fourier coefficient with respect to the 
pitch stirring frequency highly dependent on pitch stirring 

frequency. versus radial position at different excitation 

frequencies Is plotted In Figure 14b. The curves consistently peak at 
approximately r/R .8. 

The dynamic Induced flow for a pitch stirring frequency of u ■ .933 
from Figure I4d Is plotted In Figure 15. A linear radial distribution 
assumption Is not well satisfied. Straight line distributions that 
produce the same hub moment as the experimental distributions were 
calculated at all pitch stirring frequencies. The straight line In 
Figure 15 is the equivalent radial d)mamlc Induced flow for the 
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plotted experimental Induced flow. Theae values are plotted in 
Figures 16, 17, and 18 for 6^ •• 2, 5, and 8 degrees. 

Also plotted in Figures 16, 17, and 18 are the experimental 
frequency flapping responses. The theoretical curves are from the two 
parameter Identification test results described In Section 6.3. 
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Figure 14b 
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Figure 14c 
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6.2 DATA LENGTH STUDIES 

Most of our early failures in developing the parameter Identifica- 
tion methods were due to using Insufficient data lengths. Therefore 
considerable effort was devoted toward determining how much data 
should be used to Identify selected numbers of parameters for different 
models. Data length studies were made of an arbitrarily chosen test 
condition “ 5°, = 1.18, h/R - 1.02) for several different models. 

Two different data length studies were made for each mode^.. The 
first consisted of finding the shortest data length to give parameter 
convergence and then Identifying the parameters using Increasing 
amounts of data to the end of the data set. The Interval of data 
length was approximately .25 rotor revolutions. This produced plots 
of parameters Identified versus data length. This direct study was 
compared with the method of optimal data utilization developed In 
Reference 42, which was built Into the identification algorithm 
(Section 5). The whole data length was used to identify the parameters. 
On the Iteration following parameter convergence the program inverts 
H every 15 time steps (Section 5.3.3), thus the data length study 
Is a by-product of the Identification process. 

The results of the data length studies for the quasisteady wake 

model (12) are presented In Figure 19. A* does not converge to a 

—1 *s 

constant value with data length nor does (M )^ level off with 
Increasing data length. Thus no particular data length Is optimal for 
Identifying A* . The reason for this Is related to the fact that the A* 
model has basic limitations, see Figures 17 and 27. An arbitrary data 
length of 9.0ir was used In Identifying the A* quasisteady wake model. 

The results of a 2 parameter (L, t) data length study of the 
unsteady wake model, which contains 4 parameters (P2, A, L, t) are 
presented In Figure 20. P2 .3924, the structural parameter. Is 
assumed known. A =* .385 comes from blade element momentum theory. 

After t “ lOir both the square root of the diagonals of M ^ and the 
identified parameter values have reached a steady state condition. Both 
studies produce the same optimal data length. A data length of lOir 
Is used In Identifying the L and t first order unsteady wake model. 

The results of a 3 parameter (A, L, t) data length study of the 
unsteady wake model are presented in Figure 21. Although A Is 
Identified as .385, the same as the specified value of A In the 
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2 parameter modelt it takes at least 12ir for the parameter values and 
the square roots of the diagonal terms of to reach steady-state 
values. A data length, of 12ir Is used In Identifying A, h, and t. 

The results of the 4 parameter (P2, A» L, t) data length studies 
are presented in Figure 22. Neither the parameter values nor tha 
square roots of the diagonals of M~ reach steady state by the end of 
the data set, although they are apparently very close to doing so. 

The Identified value of P2 Is .38 which Is very close to Its true 
value .39. The four parameter unsteady wake model was not Identified 
except for the one collective pitch presented In Figure 22. 

Figures 20, 21, and 22 show that as the number of Identified 
parameters is Increased for a particular model, longer data lengths 
must be used to Identify the model. They also show that both data 
length studies produce approximately the same result. Since the M ^ 
study is on the order of 15 times less time consuming than the 
Identified parameter data length study, it is to be preferred over the 
Identified parameter study. 

From Figure 27 one may suspect that the reason for the improve- 
ment in the identified parameters when using data beyond « Sir is that 
the input becomes more like a constant amplitude oscillation. 

However, as can be seen from Figure 30, a short Input transient results 
in a response that extends over 8r. beyond the end of the Input. As 
discussed later identification from this short input transient using a 
data length of lOir, yielded almost the same parameters as the input 
shown In Figure 22. 
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6.3 IDSHTIFIED PABAMBTESS 

Figures 23 , 24, 25, and 26 presoit identified parameters of the 
perturbation models for dynamic induced flow from progressing pitch 
stirring tests. Plotted are the fit factors, RR (22), versus collec- 
tive pitch and the pa r a m eter values versus collective pitch for each 
data set. 

Figure 23 presents the results of Identifying the quasisteady 
wake model, (12), from Transient Tests Sets 1 and 2 of Table 4. 

Figure 24 presaats the results of identifying 2 parameters 
(L, t) of the 4 parameter (A, P2, L, t) first order unsteady wake 
model from the Transient Tests Sets 1, 2, and the first part of 3 
of Table 4. The meaoured parameter P2 varied slightly from experiment 
to experiment and was adjusted to the correct value for each data set. 
A = .385 was assumed constant, and Is the calculated value of A at 
a collective pitch of 5 degrees. 

Figure 25 presents the results of identifying 3 parameters 
(A, L, t) of the 4 parameter first order unsteady wake model from the 
Transient Tests Sets 1, 2, and part cZ 3 of Table 4. 

'Hie average fit factors (root mean square error over all the 
data set), RR, for tue three models are shown In Table 8. Also shown 

Table 8 

Averaged Fit Factors 


Model RR Conq>arison 

a' (12) .0774 2.19 
L, t(8) .0442 1.25 
A,L,t(8) .0353 1.00 
A,L,t,H(9) .0346 .98 


in Table 8 is the average fit factor of identifying 4 parameters 
(A, L, T, H) of the 5 parameter (P2, A, L, t, H) first ordqr unsteady 
wake model with wake zross-roupling (Equation (9)). The A, L,t (8) 
model is 25Z better then the L,t (8) model and 119% better than the 
A* (12) model. The A, L, H (9) model is only a rlight, 2%, 
Improvement on the A, .. t (8) model. Figure 24a shows that the fit 
factors for the A* (12) model get slightly better with increasing 
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collective pitch. Figure 25a shows that the fit factors for the L» t (8) 
model are best arouad 6^ « 5° and worsei wit**, both increasing or 
decreasing collective pitch. Figure 26a shows that the fit factors for 
the A, L, T C8) model does not vary with collective pitch. Thus the 
identified A, L, t (8) perturbation model shown in Figures 2Sa and b 
is clearly the best of the 3 perturbation models in Figures 23 to 25, 
since it has the best average fit factor and no fit factor variation 
with collective pitch. 

When the wake cross-coi^llng term, H, from Equation (9) was 
Ideitlfied along with A, L, and t, the identified values of A typically 
increased by .5%, the identified values of L typically decreased by 
1.5Z, and the idmtifled ''alues of t typically Increased hy 7Z. The 
Idoitified value of wake cross-coupling was always in the direction of 
rotor rotation and was not an obvious function of collective pitch. 

The average value of H for all data sets was .15. The arctangaat of 
.15 is 8.5 degrees. Therefore the induced flow from rotor drag in 
the direction of rotor rotation is approximately 2Z of the rotor 
velocity. 

The 5 parameter 2nd order unsteady wake model CF2, A, L, T^, T 2 ) 
of Equation (11) could not be successfully identified from the progress- 
ing pitch stirring transients. 

Figure 26 presents the results of identifying the A, L, f (8) 
model from Tremsient Tests Set 3 of Table 4. These results show 
that the Id^tified values are not functions of A, L, and t are all 

rotor aerodynamic or aerodynamic terms and should be independent of 
which is strictly a structural property of the ro'.or. 

The simulation studies of Reference 38 determined that progressing 
pitch stirring transients should be used over regressing pitch stirring 
transients. Table 9 presents comparisons of identification results, 
fit factors, and the square root of the related diagonal term of 
from different experimental transients. Regr stands for a regressing 
pitch stirring transient test, prog stands for a progressing pitch 
stirring transient test, and pspring stands for a progressing spring 
driven pitch stirring (Section 2.2 1) transient test. Two different 
regressing excitations were used for identifying parameters. These 
excitations can be seen in Figures 28, 29, and 30 of Section 6.4. The 
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_i u 

smaller the (M ) terms are» the better the eKcitatlon Is (Section 5.2). 

Table 9 shows that the progressing pitch stirring excitation is 

someirtiat better than the slow regressing pitch stirring excitation, 

~1 % 

since the QI 44 ) 'o for the progressing cases are better than those of 
the regressing cases and A is the strongest parameter. Both 
these excitations produce almost the same parameter values. The fast 
regressing pitch stirring excitation will not converge for the 3 
parameter model and does not produce good results for the 2 parameter 
model. 

Table 9 shows that the progressing spring excitation is not 
sufficient for identifying the three parameter model, but does better 
than the progressing pitch stirring excitation for identifying L and 
T if A is specified at the identified value from the progressing 
tests. The L valu.^s identified are almost equal, but the t values are 
about 20 Z less for the progressir.g spring excitation. 

Figure 23b shows how A* , the quasisteady wake parameter, should 
vary with ground effect and collective pitch according to b^^de element 
momentum theory. A* was calculated from A* » A/(li-AL) where A was 
assisned to be .385. L was calculated from (14). was calculated from 
blade element momentum theory corrected for ground height variation by 
the relaMor:. found in Reference 36. The e:q>erlmental curves for the 
smallest (h/R .78) and largest (h/R » 1.28) ground heights from 
Figure 23a are also shown. For h/R » 1.28 the eiq>erimental and theore- 
tical curves are close, however they deviate increasingly with decreasing 
ground height. The e^erlmental curves show A* increasing with 
decreasing ground height. 

The three average induced flow values measured with h/R = 1.28 
shown in Table 7 (Section 6.1) can be used to calculate L and t from 
(14). The results of these calculations are shown in Table 7. By 
comparing these values with those of Figure 24b for h/R ** 1.28 it is 
seen that there is reasonable agreement between the theoretical values 
of L and t and those Identified from transients. 

Figure 25c presents theoretical values of L and t plotted as a 
function of ground height and collective pitch. L and t were calculated 
from (14). The average induced flow values were calculated from blade 
element momentum theory corrected for ground height variation by the 
relations found in Reference 36 (Section 6.1). These curves show 
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excellent agreement with the experimental curves of 25b for L above 
collective pitch settings of 5°. At low collective pitch settings 
the identified L values are lower than the theoretical ones. The 
identified t values agree with the theoretical ones at very low collec- 
tive pitch. As collective pitch increases the identified t values are 
proportionally larger than the theoretical ones. The theoretical and 
identified values both show the same effect of ground height on L and 
T. As collective pitch increases, disk loading Increases, and tip 
losses (B) should increase. This would cause A to decrease. Figure 
25a shows the identified A increasing until the collective pitch 
reaches about 5 °. Above 5 ° A decreases as predicted by the increasing 
tip loss factor. 


Table 9 

Transient Identification Comparisons 


A, L, T (8) model 

A 


T 

TX 

RR 

oil =1.18 






h/R =1.28 

regr slow 
psprlng 

.356 .0061 
.277 .0092 

4.72 .10 
3.43 .10 

7.80 .35 
7.22 .54 

.0372 

.0367 

=1.02 

prog 

regr 

.356 .0040 
.381 .0064 

4.66 .13 
5.50 .10 

8.50 .40 
7.80 .29 

.0400 

.0357 

= .78 

prog 

regr 

.359 .0040 
.358 .0056 

4.71 .14 
5.58 .10 

9.25 .46 
8.31 .42 

.0468 

.0343 

u)2=1.26 






b/R =1.28 

prog 

pspring 

.355 .0040 
.455 .0338 

4.32 .14 
4.94 .16 

7.86 ,46 
5.28 .42 

.0354 

.0616 


L, T (8) model 
<oj“l. 18 


h/R =1.28 regr slow 

.360 

4.73 .10 

7.69 .32 

.0372 

regr fast 

- 

11.15 .90 

19.76 2.25 

.0538 

pspring 

- 

4.02 .06 

5.30 .24 

.0382 

wi=1.26 





h/R =1.28 prog 

- 

4.57 .15 

7.71 .40 

.0348 

pspring 


4.59 .14 

5.98 .47 

.0623 
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Fxgure 25a 

Fif FaCkOr and A Versus 0^ for the Three Parameter Unsteady Wake tiodel 
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Figure 25b 

L and t Ver^ius 9 for the Three Parameter Unsteady Wake Model 
o 
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Figure 25c 

Theoretical Values ot L and t Versus 8 

o 
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6.4 PREDICTIONS 

Figures 16, 17, and 18 in Section 6.1 present frequency response 

test prediction results. The crosses are the measured data points 

discussed 1.. Section 6.1. The solid lines in Figures 16, 17, and 18 

are predictions for the L, t (8) model of Figure 24. L and t values 

for 0^ » 2, 5, and 8 degrees and h/R « 1.28 were picked from the solid 

line curves in Figure 24. A was assumed to be .385. The dashed line 

in Figure 17 is a prediction for A « .385 where the dynamic rotor inflow 

is omitted. The dot-dashed line in Figure 17 is a prediction for the 

A* (12) model of Figure 23. The A* value for 6 » 5*^ and h/R ■* 1.28 

o 

was picked from the solid line curve in Figure 23* 

Omitting the dynamic inflow produces poor predictions of the 
experimental frequency flapping responses (Figure 17)# The quasisteady 
wake model produces better prediction results (Figure 17) but cannot 
predict the measurement trends as a function of frequency* Though no 
flow measurements are used to establish the solid curves (the L, t (8) 
identified model) in Figures 16 > 17, and 18, and though the presented 
flow data (crosses) are derived from inaccurate measurements (see 
Figures 14c and d) the agreement in the dynamic flow Fourier coeffi- 
cients is reasonable. The experimental frequency flapping responses 
agree even better with the L, t (8) identified model predictions. 

Figures 27 through 30 present real time prediction results. 

The dots are the measurements. Every 5th measurement is plotted. The 
experimental input for each transient excitation is plotted above the 
experimental and predicted responses. The root mean square fit 
between prediction and measured response over the length of the data 
is noted on the figures. Each of the experimental responses was used 
for identification purposes, consequently the data bias terms (Ej and 
'^II) are known. Each prediction includes these bias terms. The 
parameter values for each prediction are determined by picking them 
from the curves (Figures 23, 24, and 25) for the model being ^ered. 

Figure 27 presents prediction results, using the A, L, t 
model and the a' (12) model, from a progressing pitch stirring transient 
input. This is the transient input used to produce the models shown in 
Figures 23, 24, and 25. The predicted dynamic downwash (from the 
A, L, 1 (8) model) is plotted along with the flapping responses 3 nd 
the experimental input versus rotor azimuth angle. The quasisteady 
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model does not predict the measured response very well. The first 
order unsteady wake model predicts the measured response adequately. 

Figure 28 presents prediction results, using the A, L, t (8) 
model, from a slow regressing pitch stirring excitation. The predic- 
tion is excellent. 

Figure 29 presents prediction results, using the A, L, t (8) 
model, from a fast regressing pitch stirring transient input. 

Figure 30 presents prediction results, using the A, L, t (8) 
model, from a short progressing pitch stirring excitation. This is 
the excitation produced by the spring loaded drive described in Section 
2.2.1. There is a two per rotor revolution component in the measured 
response. Aside from this, the prediction is excellent. It has also 
been shown (Table 9) that the L, t (8) model can be identified from 
this transient. 



1.6^ e = 1.5°, h/R a 1.28 
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Predlcced with 
A, L, T (8) model 
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7 , CONCLUSIONS 

Since the objective was to gain a better understanding of rotor 
aeroelastlc phenomena from dynamic rotor tests using both transient 
tests O^th state variable and parameter Identification evaluation 
techniques) and frequency response tests (with Fourier analysis evalu^ 
at Ion techniques), there are two types of conclusions; the first type 
regarding methods used; the second type regarding the results for the 
selected aeroelastlc problem* 

7.1 CONCLUSIONS REGAKDING METHODS 

1.1 Transient tests evaluated by a simplified version of the maximum 
likelihood state variable and parameter Identification method 
reduced model testing time by a factor of 80 over harmonic tests 
evaluated by frequency response methods for each rotor condition 
studied. 

1.2 The parameter Identification method provided the model parameter 
values directly while model parameter values must be determined from 
curve matching from frequency response tests. 

1.3 The method of optimal data utilization based on the time history 
of the inverse of the sensitivity matrix can be effectively used to 
prevent inaccuracies from Insufficient data utilization* 

7.2 CONCLUSIONS REGARDING ROTOR DYNAMIC INFLOW EFFECTS 

2.1 Transient cyclic pitch stirring Is a simple and effective test 
technique to extract linear perturbation models Including dynamic 
Inflow effects from blade flapping measurements* The perturbation 
model used Is only valid In the frequency regime typical for flight 
mechanical or ground and air resonance phenomena* 

2.2 The analysis using identified rotor dynamic Inflow parameters 
Is In good agreement with transient and frequency response test 
results not used In the Identification; thus the extracted perturba- 
tion is validated within the stipulated frequency range* 

2.3 The extracted perturbation model is in reasonable agreement with 
a simple dynamic inflow theory especially above collective pitch 
settings of 5 degrees, thus confirming the previously deduced 
phenomenon of a substantial reduction of the regressing flapping 
mode damping caused by dynamic rotor inflow. 

2.4 The ground proximity within the tested range (h/R • • 78 to 1.28) 
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has a pronounced effect on the dynamic rotor downwasb parameters 
which become larger with decreasing ground-rotor distance. 
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9. APPENDIX 


9.1 NOMENCLATURE 


lu) 


A 

A* 

A" 

B 

B 

F 

H 

H* 

I 

J 

K 

L 

L* 

L' 

M 


a b^y/ 8 aerodynamic blade conscant, (4)* 

« A/d-t-AL) aerodynamic blade constant assuming 
quasisteady vake» (12) 

aerod}mamlc wake constant » (4) 

blade tip loss factor 

measurement equation error covariance » (19) 
first Fourier coefficient of frequency at 
wake cross coupling term, (9) 

- H/t, (9) 

>■ apcR^/y blade moment of inertia 
scalar cost criterion, (17) 

1024 words 

dynamic Induced flow gain, (5), (13) 

« AL/t, (8) 

- AL/tJ, (11) 

Information or sens:.tlvity matrix, (20) 


%wi, ”aA’ 


«b 

N 

P2 

R 

RR 

T 

T, 


diagonal terms of M 
blade moments 

number of measurement vectors 


> (u^ - 1) non-dimensional non-rotating blade flapping 
natural frequency (8) 

rotor radius 

fit factor, (28) 

- T2/T^, (11) 

acceleration term of the second order unsteady wake 
model, (10) 


* Numbers in brackets here and in the text refer to equations 
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T' 

W 

m 

a 

b 

c 

fOl 

h. 

hC-) 

n 


rate cent of the second order unsteady wake model, (10) 

a + AL)/tJ, (11) 

hot wire anemometer wake measurement 
average wake measur^nent 
blade section lift slope 
number of blades on the rotor 
blade cord 

function of variables In parentheses In the system 
equation, 0-5} 

rotor to ground plate distance 

function of variables In parentheses In the measurement 
equation, (16) 

rotor speed 


r 

u 

V 

w 

X 


y 


^Im’ ®IIm 

g 


\l>, t 



X 


k 


distance from rotor center along blade span 

input control vector, (IS) 

measurement noise vector, (16) 

system or process noise vector, 05) 

state vector, 05) 

measurement vector, (16) 

flapping angle of the kth blade, positive up 

multiblade flapping coordinates: longitudinal and 

lateral cyclic flapping 

measurement bias values 

multiblade flapping coordinate measurements 

estimated values of flapping angles 

rotor azimuth angle 

azimuth angle of the kth blade, (5) 

inner shaft or cyclic pitch azimuth angle, (1) 

a constant value, (2) 

generalized vake velocity, (4) 
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(II 


“a 



a 

T 


cyclic pitch stirring excitation frequency in the 
rotating system 

cyclic pitch stirring excitation frequency in the 
space-fixed reference system 

non-dimensional rotating blade fleqiping natural frequency 

innovation vector et time j, (18) 

V uniform, longitudinal and lateral perturbation induced 
Inflow coniKmeota 

r weighted averse Inflow 

air density 

blade Lock nund>er 

vector of unknown parameters, (IS) 

pitch angle of the fcth blade positive when leading edge 
is up 

^1* ^11 cyclic and left cyclic pitch angles, 

respectively, (6) 

rotor solidity ratio (b*c/irl0 

dynamic Induced flow time constant, (5) 


T* = a + AL)/t, (8) 


Superscripts 

• time derivative 

estimate 

T 

transposed matrix 
~ bias value 


Subscripts 

i 

k 

m 

P t 
t 


value for Jth data sanple 
value for kth blade 
measured variable 
pretranslent 
transient 
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o, I. II 

lo) 


multlblade variables 
first Fourier coefficient 
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OS/0 


0O0S 

0003 

0004 

0005 

0006 
0007 
0010 
0011 
0012 


0013 


0014 


0015 

0016 

0017 

0020 


0021 

002? 

0023 


00?4 

00?0 
0027 
H030 
W0 31 
00 3? 


FORTRAN IV 3,03 APR 24 J 077 PARF ONF 

C HVPRT 4 IS A PARAMETER IDENTIFICATION 
C PROGRAM DESIGNED TO IDENTIFY 4 TO 6 
C PARAMETERS (L* TAU^SlASl ,BI A$ 2 , A,P 2 ) 

C FORM MOVERING TRANSIENT TESTING, 

niHENSION BM( 3 | 500 ),RJ( 6 ) ,Y( 42 ),RM(S, 6 )fRRR( 2 | 2 ) 
DIMENSION PRMT( 3 ),Gt 2 ,l 5 ),PR( 20 ),Z( 6 , 2 , 15 ),RR( 6 ),AU 7 ( 42 ) 
COMMON M,MM,V,AiJX,PP#NP,NSD, 0 * 1 ,G,Z,RM,RJ,RR,RRR 
external FC 4 THV, 01 ! 4 THV,mINV 

c read in data 

REAO( 5 , 20 )IJK,PNT 8 ,F 1 RZO,PITO 

20 F 0 RMAT(I 4 , 3 FtS,R) 

DO 10 I«1,S 

READ ( 5 , 21 ) (MM (I, J), Jr 1 , UK) 

21 FORMAT( 10 F 7 , 3 ) 

C rjKvNUMBER OF MEASUREMENT VECTORS 
C PNTS«NUMPER OF MEASUREMENTS PEP ROTOR REVOLUTION 
C FIRZOsNUMBER OF MEASUREMENTS INTO THE DATA SET 
C THAT THE FIRST ZERO AZIMUTH ANGLE OCCURS 

C PITObINITIAL condition of THE PITCH AZIMUTH ANGLE 
10 CONTINUE 

C 3 M( 1 , )*PITCH AZIMUTH 
C 5 M( 2 , )« 8 ETAl 
C RM( 3 , )»PETA 2 
C ******************** 

c NpBurF parameters 

NP*5 

c ******************** 

C 4 YEILDS L»TA(',MIAS 

C 5 YEILDS L»TAI.',PIAS, A 

C 6 YEILDS LiTAU.eiAS, A,P** 2 ,«l, 

NS 0«6 

C nSO*DImEnSION df state eon, 

NIC«(NP-n*NSD 

C NITbPOF initial rONOITlO^S 
PRMTf DtFIRZO 
PRMT ( 2 ) ■ 6 , 8318531 /PNTS 

C BRMT( 2 ) IS THE TIME UNIT BETWEEN SAM; UE POINTS 
C IN RADIANS AND THE INTERVAL OF INTEGRATION, 

RR( 1 )b 400 , 

RR( 2 )» 0 , 

RP( 4 )» 400 , 

C RRfl- 4 ) IS THE COST FliNCTION NORMALIZATION MATRIX, 

C ******************** 

C INITIAL parameter ESTIMATES ANP/QR SPECIFIED 
C PARAMETER VALUES, 

A 8 , 40 P 

XL« 4 , 

TAU» 4 , 

P ?*,3924 

B ««,976 

MRsPR**4, 

XKAPPA«,1 

r ******************** 

r PRfl- 20 ) IS THE parameter ASSOCIATEO VECTOR 
C PASSED TO SUBROUTINES, 
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08/8 FORTRAN IV 8.03 


APR 24 197/ page two 


0033 

0034 
0039 

0039 
0037 

0040 


0041 

0042 

0043 

0044 

0049 

0049 

0047 

0090 

0091 

0092 


0093 

0094 


0099 


0056 

0097 

0060 

0061 

0062 

0063 

0064 

0069 
0066 
0067 

0070 

0071 


0072 

0073 

0074 


PR(6)«P2 

C P2«P SQUARED NINtlS 1, 

PR(1)«XL 

PR(2)«TAU 

PR(9)«A 

PR(t8)«8XN(PlT0) 

PR(in«COS(PXTO) 

c ************* f ft 

C INITIAL BIAS VALUES 

PR(3)».0 

PP(4)««0 

C ******************** 

HRITE(3«22)TJK 

WRITE(3f03)A>XL,TAU,P2 

22 FORMATC19X, *6 OF 0EA8UREHENT8«» , 18) 

23 FORMATOX, ’INITIAL ESTIMATES’ »3X, *A«SF7«3,3X, 

’ ’L«’|P7,3,3X,’TAU*’,F7.3,» P2*’,F7.4) 

WRITE(3,l7nPR(3),PR(4) 

171 FORHAT(20X,«Bt 6IAS«’ ,F8.9, ’ B2 R1 ASp’«F8,9) 
WRITE(3,172)BB,XKAPPA 

172 FnRMAT(20X,*B**4,«’,F8.5,’ K APPA»’ ,F8,9) 

C ******************** 

C CALCULTAE MULTIPLVING FACTORS 

n IKK*15«THE NUMBER OF DATA POINTS USED IN IDENTIFICATION. 
IKKb20 

C ******************** 

24 F0RMAT(T41 

C ******************** 

C SPECIFIES THE ROF PARAMETER ITERATIONS. 

C ITERATION LOOP 

on 60 MM»l,9 

c ******************** 

C SET STATE VECTOR INITIAL CONDITIONS TO ZERO, 

00 70 I«1,NIC 

70 Y(I)«0, 

C ZERO THE PARTIAL J/PARTIAL PARM AND THE M ARRAYS, 

00 71 I«1,NP 
RJ(I)«0, 

00 71 J«1,NP 
RM(I,J)«0, 

71 CONTINUE 

C AVOAOR IS THE AVERAGE OF THE DIAGONAL OF B 
AV0AGR«(RR(l)^RR(4))/2. 

RRRrunaAVDAGR 

RRR(2,2)BAVnA6R 

RRR(1,2)"0, 

RRRf2,n*0, 

r «WR(1,1)«RR(1) FOP normalizing COST WITH B 

C R«R(2,2)bRR(4) for NORMALIZING COST WITH 6 

C RRP(1 ,2)*RR(2) FOR NORMALIZING COST WITH B 

C WRP(2,1)»RR(2) " " nun 

on 72 I«l,4 

RR(n«0, 

72 CONTINUE 

c COMPMTF nEw parameter ASSOCIATED VECTPR FOR EACH 
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08 /B 


0075 

0075 

0077 

0t00 

0101 

0102 

0103 

0104 

0105 

0106 
0107 
0110 

0111 

0112 
"f • ^ 

0U4 

0111 

0110 

0117 


rORTRAN IV 3,03 APR 24 1177 PAfiE 

C PARAHETEw ITERATTOM, 

PR(20)«PR(n/PR(2) 

PR(l9)»fPR(5)*PR(l)+l,)/PR(2J 

PR(10)«PR(9)*PR(1)/PR(2) 

PRC17)«1,/PRC2) 

PR(I6)«PRf5)7PR(2) 

PR(151«PR(9W(PR(5)*4,*,028*XKAPPA)/6B 

PR(14)»PR(15)*PR(1)/PR(2) 

C SET UP working vector POR INTEGRATING SUBROUTINE. 

DO 73 I«1,NIC 
73 Auxrn«0. 

C OUTER integration UOOP 
00 00 N«1,IKK 
NIJHSTP*15 
N0IM«NIC 

C INTEGRATE THROUGH 15 TINE STEPS, 

CALL RK4THVCPRPT,NIJMSTP,NPIM,FC4THV,0U4TMV) 

C RETAIN PARTIAL SUMS OF INTEGRATION RESULTS 
CALL SEN4TS(IKK,MTNV) 

60 CONTINUE 

CONTINUE 
80 CONTINUE 

CALL EXIT 
Er40 


TMRFE 
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os/e 


0flfl2 

0003 

0004 


0005 

0006 


0007 

0010 

0011 

0012 

0013 

0014 

0015 

0016 
0017 
0020 


0021 


0022 

0023 

0024 

0025 

0026 

0027 

0030 

0031 


0032 

0033 


FORTRAN IV 3,03 APR 24 1<)77 PAGE ONE 

C fC4TNV IS THE FUNCTION SUBROUTINE TO BE USEO IN 
C IN THE EVALUATION OF H0VFRIN6 PARANETER 
C IDENTIFICATION. IT IS CALLED BY RK4THV. 

SUBROUTINE FC4THVrN,DERVf J) 

DIMENSION Y(42)fDlRV(42)*BN(3«500)fPR(20),AUX(42) 
COMMON M,HN.Y*AIIX*PR,NP,N 8 D,BN 
C I 18 THE TIME AT THE St ART OF THE INTERVAL OF 
C INTEGRATION. Ilal^l IS THE TIME OF THE ENO OF THE 
C INTERVAL OF INTEGRATtOl. 

I«(M*n*lS«N 

II«(M«1)*1S4N«1 

C COMPUTE LEFT AND FOR4AR0 CYCLIC PITCH INPUT 
C FROM THE PITCH AZIMUTH ANGLE* 

C SUBTRACT THE INITIAL CYCLTC PITCH INPUT ANGLES 
C (PR(t 0 ,Un* 

60 TO (10,tliU«l21tJ 
12 THAl«.C08(BM(l,tIl)fPR(ll) 

THA2P»8IN(BM(t»tI1)4PR(10) 

GO TO 14 

11 THAla((«Cn8(BMM.|n«C08(BM(l,inn/2.)«PR(in 

THA 2 «(*SlN(RHdn 1 ) 4 StN(BM(l,in))/ 2 *aPR( 10 ) 

GO TO 14 

10 THAia»Cn 8 (BM(l,in«PR(ll) 

THA2»«SIN(BM(ttI))«PR(10) 

14 CONTINUE 

C >illa THE number OF DIFFERENTIAL EQUATIONS TO 

C BE SOLVED. 

NTIaN80*(NP«n 
C PR(20)aL/TAU 
C PRdOJalAaL^D/TAU 
C PR(lB)aA*L/TAU 
C PR(17)al./TAU 
C PR(16)aA/TAU 
C PR(15)«A« 

C PR(14)aA»*L/TAU 
C PR(ll)a*THETAt AT Tap 
C PR(l0)a « 2 " " 

C «*••****«*** 

C 00 LOOP THAT SOLVES THE KERNAL OF THE SYSTEM 

c equations and the sensitivity equations. 

00 IS I«lfNII*NSD 
0ERV(T)aV(Itl) 

nFRY(I + n«*PR(6)*V(I)«PR(l5)*V(l4n«Pfin5)*Y(Ia2) 

' ^2.*V(I^3UPR(5)*V(I44) 

nERV(T+2)aVfIt31 

nERY(lY3)tPB(l51*Y(n^2.*Y(I + n-PR(6)*Y(I^2) 

• -PR(151*V(I*3)*PR(5)*Y(I45) 

DERV(I»4)8PN(r41*V(l4ntPR(14)*V(Ia2)-PR(l9)#Y(Ia4) 

OFRY(It5)a«PR(14)*Ym*PR(14)*Y(I*3)*PRfl0)*Y(If5) 

15 continue 

c ********«• 

C AOn ON THE FORCING FUNCTIONS TO THE DIFFERENTIAL 
c equations 

0FRY(2)«PFRY(2)^P»(5)*THA2 

nFRY(4)8DERY(4)*PR(5)*THAl 
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OS/4 

0P34 

B035 

0036 

0037 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0050 


0051 

0052 

0053 

0054 


FfJOTRAN IV 3,03 APP 24 107/ PAfif TWO 

ocoy (5) aoeov (S) * pp ( i ai • th a 2 
OERY (6) «DERy ( 6) aPR O 0) *THA1 
C (2*6) ARE SYSTEM EQUATION FORCING FUr'CTIOMS, 

OERy(ll)»OCPY(ll)aPR(16)*Y(2)fPR(16)*Y(3)*PR(10)*Y(S) 

* •PR(16)*THA2 

OERY(l2)0OERY(l2}«PPn6)*Y(napRri6)*Y(4)«PRn0)*V(6) 
’ ♦PR(16)*THA1 

0ERY(17)f0ERY(17J-PRfl7)*0ERY(5) 

OERY(lB)0OERy(10)*PRCl7)*OERY(6) 

XP(NP.5)20*21»22 
22 0ERY(26)»nERY(26)*Y(l) 

0ERY(2S)«0ERY(26).Y(S) 

21 OERY (20) aOERY (20) «Y (2) .Y (8) ♦ V (5) 40 

nERY(22)*nERY(22)4Yri)«Y(4)4Y(6)»THAl 
OERV(23)«f>ERY(23)4PR(20)*Y(2)4PR(20)*Y(S)«PR(20)*Y(5) 

* •PR(20)*TMA2 

OERY(24)*nERY(24)«PR(2e)*y(l)4PR(2O)*V(4)«PR(20)*Y(0) 

* «PR(20)*TH41 

(n»24) ARE SENSITIVITY EQUATION FORCING FUNCTIONS. 

20 CONTINUE 

RETURN 
STOP 
Eno 
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0002 

0003 

0004 
0009 


0006 

0007 


0010 

0011 

0012 

0013 

0014 

0019 
0010 
0017 

0020 
0021 
0022 

0023 

0024 

0029 


0026 

0027 

0030 

0031 

0032 


0033 

0034 


0039 
0030 

0037 

0040 

0041 

00 4 ? 


C 0U4THV IS THE OUTPUT SUBROUTINE CALLED PROM 
C RK4THV TO BE USED WITH HOVERING DATA ANALYSIS, 

SUBROUTINE 0U4THV (lI,DERY,NOIM,PRNn 

DIHENSION 0ERY(42),Y(42),PRNT(S),0(2,19),Z(0,2,19) 

DIMENSION BH(9,9001 ,PR(20) * AUX(42) 

COHMON M,HH,Y,AUX|PR,NP,NS0,BH,6,{ 

C BEST ESTIMATES OP LEPT AND FORWARD TILTS AT THE END 
C OF THE INTERVAL OF INTEGRATION ARE SAVED, 

6(1«II)«Y(1) 

G(2,XI)tY(3) 

C ********** 

C PARTIAL J/PARTIAL PARAMETER VALUES ARE SAVED, 
Z(l,l,Xl)*Yf7) 

Z(lr2»lI)"Y(9) 

Z(2,lf in«Y(13) 

Z(2,2,II)«Y(19) 

Z(3,l,II)il, 

Z(3,2,2I)«0. 

Z(4,l»in*R, 

Z(4,2,II)*1, 

XF(NP«9130,31,3? 

32 Zr6,l,II)*Y(29) 

7(6,2, II)*Yf27) 

3t Z(5,l,in*V(19) 

Z(9,2,Zt)*Y(2t) 

C ********** 

30 CONTINUE 
C ******************** 

C INPUT, OUTPUT, ANO STATE VECTORS AS A 
C PUNCTION OF TIME CAN BE GENERATED DURING 
C ANY ITERATION BY SUBSTITUTING 25 FOR 40, 

GO TO (40,40,40,40,40,40)|HH 
C ******************** 

40 RETURN 

25 l*(M«l)*15*II*l 

THAl*«CnS(BM(l,n)*PR(ll) 

THA2*i.SIN(BMn,n)*pR(l0) 

C Y| AND Y3 ARE THE JEST ESTIMATES OF THE 
c tilting states, 

Yl*Y(l)*PR(3) 

Y3*V(3)*PR(4) 

C SIM and BXN2 ARE THE PfiST ESTIMATES OF THE 
C JAKE STATES, 

BIN1*V(9) 

BIN9*Ytfl5 

C SMf 2 ,I) ANO RMf3,n ARE THF TILTING MEASUREMENTS, 

WRITE(3,24)BM(2,n,Vl,8lNl,BM(3,n,V3,BlN2,THAl,THA2 
24 F0RMAT(4)(,flF8,3) 

RETURN 

End 
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0«02 

0003 

0004 

0005 

0006 

0007 


0C10 

0011 

0012 

0013 

0014 

0015 

0016 
0017 

v’P20 

v \ v >2\ 

0P22 

0023 

0024 

0025 

0026 
0027 

0030 


0031 

00.32 

0033 

0034 
003?^ 

0036 

0037 

0040 

0041 

0042 

0043 


0044 

0046 

0046 

ifiCA ? 

005'-’ 


FORTRAN IV 3,33 APR 24 1**77 PaGF 0* t 

C 8F» 4TS COMPUTES THE PAPTAU SUMS OF The 3F. *STT1VTTV 
C 3EOUFNCES, 

SUBROUTINE SE'UTSf I»<K,MINV) 

DIMENSION OH(4),Gr2.15),Z(6,2,l5),R^(6,6),eM(3,S06),RRR(2f2) 
DIMENSION RI(6),LV(6),MV(6),PR(20)«Y(42),6RM(36)fRR(6}«AllX(42) 
dimension 0TAGH(6) 

COMMON H,MM,Y,AUX,PP,NP,NS0,eM,G,2,RM,RJ,RR,RRR 

II«(M«n*15 

c FORM INOVATION SEQUENCE FOR THIS 15 SAMPLES 6, 

c ALSO calculate b matrix, 
no 12 J«1,16 
III#II4J*1 

C COFPUTE INOVATIOM seoupnef 

Gn,j)«6Mf2,iin-r.(i,.n*»PR(3) 

G(2,JI«HM(3,Iin»0(?,Jl-PR(4) 

c compute b matrix 

RR(l)*RRn)+G(l,JI*Rfl,J) 

RR(?)«RR(a)»G(l,J)*R(2, J) 

RR(4)«RP(41+G(2,J)*r(2,J) 

12 CONTINUE 

C COMPUTE PARTIAL PARTIAL PRAMS, RJ, 
on 14 JJ«l,up 
no 14 J»1,1S 

R.I(JJl«RJfJJ)»6(l,J)*(Z(JJ»l»J)*RRR(l,n4T(JJ,2,JI*RRR(2»in 
» ♦G(2,Jl*(ZrjJ,t,J)*RRR(l,2)^RRR(2,2)*Z(JJ,2,J)) 

14 CONTINUE 

C COMPUTE M 

no 47 T»1,NP 

nn 47 J«i,NP 
on 47 K«l,15 

RM(I, J)«RM(I, l) + (Z(I,l,K)*RRRn,n+Z(I»2,K)*RRR(2,in*Z(J,l,K) 
' ♦(zn»l.'<l*RRR(l,2I+Z(t,2,K)*RRR(2,2n*Z(J»2,K) 

47 CONTIauE 

C SUPSTTTUTE 31 FOR 33 FOR PRINTING OUT TIME HISTORY 

C OR M inverse 

GO TO (33,33,33,33,33) . MM 
33 IF(M,F0,IKK) GO TO 31 

Gn TO 32 

C CHANGE RM TO GENFRAL MATRIX GRM, 

31 nr«p 

nn 49 J*l,NP 
no 49 I«1,NP 
UP«M8*1 

GRM(HB)*R‘<(T, J) 

49 CONTlwiJE 

NP2«NP*MP 

[. A******************* 

23 FnRMAT(4X,5F10,6) 

C ««********•*****•••• 

C TNVEPT M 

CALI Ml*vfGPM,NP,nT,LV,MV) • 

WRITE(3,24)m,OT 

24 FORmaT(5X, I3,3X, 'INVFPSF OF M^MATRIX',* 0Tm« I ,E 14,H/) 
>lRITFf3,23) (GRM(!) ,Isi ,NP2) 

•.RITE(3,26) 
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0051 

26 

F0RHAT(5X, 'SO RT OF OIAR M') 

0052 


DO 27 I»1,NP 

0093 


lLM"(I«n*NP4l 

0094 

27 

DIA6M(n«(AB8(GRH(ILM)))**,5 

0055 

0056 


WRITE (3,23) fniA6H(n,I»l,NP) 
WRITE h«29) (RJ(ntX«l,NP] 

0097 

25 

FORMAT (37, 'PART, ( J/THETA) ■ ' ,6F9,0/) 

0060 

32 

continue 

0061 


IF(M,E0,1KK) GO TO SO 

0062 


RETURN 


C 1 

CHANGE 6RM(K3 TO HATRIK RM(I,J) 

0063 

30 

NBS0 

0064 


on 50 J«1,NP 

0065 


00 50 I«1,NP 

0066 


NRaNB^l 

0067 


RH(t, J)aGRH(NH) 

0070 

50 

CONTINUE 

0071 

Q0 

FORMAT(4F10,4) 


C 

UPDATE UPDATES PARAMETER ESTIMATES AFTER EACH 


c 

ITERATION, 

0072 


on 95 I«l,NP 

0073 


OH(n*0, 

0074 


on 59 J«1,NP 

0075 


nHm«RM(t, j)*Rj(j)*0H(n 

0076 

55 

CONTINUE 

0077 


no 56 I«1,NP 

0100 


PR(n«PR(I)*OH(I) 

0101 

56 

CONTINUE 

0102 


RR(3)«RR(2) 

0103 


XN«19*IKK 

0104 


DO 34 I«1,4 

01^19 

34 

RR(T)»RR(I)/XN 

0106 


RRP»((RR(n*RR(4n/2.1**,5 


C 

7RP TS THE FIT FACTOR 

0107 


WRITE(3,6nRRP 

0110 

61 

FORMAT(/»RRP«',F20.6,/'RRii') 

0111 
01 12 

52 

WRITE(3.62) rRR(I),I«l,4) 
PnRMAT(2F20,8) 

0113 


CALL MINV(RR,2,0T,LV,MV) 

0114 

53 

FORMAT(4Fl0,71 

0115 


A«PR(3) 

0116 


TAU«PR(2J 

0117 


XL»PRCn 

0120 


P2»PRf6) 

0121 


wRITE(3,20)MM,4,XL,TAU,P2 

0122 

20 

FORNATC ITERATIONS*', II, ' A*',F8,3,* L*',F6.3,' 

' F6.3,' P2*',F8,51 

0123 


MRITE(3,2nPR(33 ,PRf4) 

0124 

21 

FORMAT(20X, »Bt 8 1 AS* ' , F0 , 9 , ' 82 BI AS* ' , F8,9) 

0125 


rftmrn 

0126 


End 


82 BIASa',F8,9) 



